source conda/bin/activate privratsky rmarkdown::render(‘./9_Subcluster_Macrophages/9_Subcluster_Macrophages.Rmd’)

Changes in myeloid and kidney cells after CLP - Analysis of 2 x 10X scRNA-seq samples from 2 pools of WT mice (3 Sham + 3 CLP): comparison of gene expression in different cell populations

Upregulation of anti-inflammatory mediators in macrophage population (SOCS3, IL-1R2, IL1rn) and downregulation of pro-inflammatory genes

indir <- "./processedData/7_DEA_C24_vs_C0"
outdir <- "./processedData/9_Subcluster_Macrophages"
dir.create(outdir, recursive = T)
library(Seurat)
library(ggplot2)
annotated <- readRDS(paste0(indir, "/3.annotated.rds"))
annotated
## An object of class Seurat 
## 24399 features across 18055 samples within 2 assays 
## Active assay: RNA (22399 features, 0 variable features)
##  1 other assay present: integrated
##  2 dimensional reductions calculated: pca, umap

Subset to Macrophages and sub-cluster

annotated
## An object of class Seurat 
## 24399 features across 18055 samples within 2 assays 
## Active assay: RNA (22399 features, 0 variable features)
##  1 other assay present: integrated
##  2 dimensional reductions calculated: pca, umap
Idents(annotated) <- "annotation.1"
macrophages <- subset(annotated, idents = "Macro")
macrophages
## An object of class Seurat 
## 24399 features across 588 samples within 2 assays 
## Active assay: RNA (22399 features, 0 variable features)
##  1 other assay present: integrated
##  2 dimensional reductions calculated: pca, umap
table(macrophages$annotation.1)
## 
##          Endo          Podo         PT-S1         PT-S2         PT-S3 
##             0             0             0             0             0 
##           LOH Distal Conv T   Conn Tubule         CD-PC         CD-IC 
##             0             0             0             0             0 
##      CD-Trans         CD-IM           Fib         Macro           PMN 
##             0             0             0           588             0 
##       B lymph            NK 
##             0             0
DefaultAssay(macrophages) <- "RNA"
macrophages.list <- SplitObject(macrophages, split.by = "sample.id")
macrophages.list
## $C0
## An object of class Seurat 
## 24399 features across 267 samples within 2 assays 
## Active assay: RNA (22399 features, 0 variable features)
##  1 other assay present: integrated
##  2 dimensional reductions calculated: pca, umap
## 
## $C24
## An object of class Seurat 
## 24399 features across 321 samples within 2 assays 
## Active assay: RNA (22399 features, 0 variable features)
##  1 other assay present: integrated
##  2 dimensional reductions calculated: pca, umap
for (i in 1:length(macrophages.list)) {
    macrophages.list[[i]] <- FindVariableFeatures(macrophages.list[[i]], 
        selection.method = "vst", nfeatures = 2000, verbose = FALSE)
}
k.filter <- min(200, min(sapply(macrophages.list, ncol)))
k.filter
## [1] 200
macrophages.anchors <- FindIntegrationAnchors(object.list = macrophages.list, 
    dims = 1:30, k.filter = k.filter)
macrophages.integrated <- IntegrateData(anchorset = macrophages.anchors, 
    dims = 1:30)
DefaultAssay(macrophages.integrated) <- "integrated"
macrophages.integrated <- ScaleData(macrophages.integrated, verbose = T)
dim(macrophages.integrated)
## [1] 2000  588
macrophages.integrated <- RunPCA(macrophages.integrated, npcs = 50, 
    verbose = T)
system.time(macrophages.integrated <- JackStraw(macrophages.integrated, 
    num.replicate = 100, dims = 50))
##    user  system elapsed 
##  74.031   0.002  74.104
macrophages.integrated <- ScoreJackStraw(macrophages.integrated, 
    dims = 1:50)
j <- JackStrawPlot(macrophages.integrated, dims = 1:50)
pdf(paste0(outdir, "/1_JackStrawPlot_macrophages.pdf"))
j
dev.off()
## png 
##   2
macrophages.integrated <- FindNeighbors(macrophages.integrated, 
    dims = 1:50)
macrophages.integrated <- RunUMAP(macrophages.integrated, dims = 1:50)
for (i in seq(0, 2, 0.1)) {
    macrophages.integrated <- FindClusters(macrophages.integrated, 
        resolution = i, verbose = F)
    macrophages.integrated[[paste0("macrophages_subclusters_res.", 
        i)]] <- macrophages.integrated$seurat_clusters
}
head(macrophages.integrated[[]])
##                      orig.ident nCount_RNA nFeature_RNA percent.mito sample.id
## AAAGGGCAGTCACGAG--C0         C0       1233          661     8.191403        C0
## AAAGTCCCAGGTTCAT--C0         C0       3888         1562     9.027778        C0
## AAAGTCCTCCGAGGCT--C0         C0       4237         1508     2.336559        C0
## AACAAAGCACCTGCGA--C0         C0       3528         1579     4.365079        C0
## AACCTGACATGCGGTC--C0         C0        884          420    19.230769        C0
## AACGTCAAGGTTCATC--C0         C0       4528         1895     5.631625        C0
##                      integrated_snn_res.0 seurat_clusters
## AAAGGGCAGTCACGAG--C0                    0               6
## AAAGTCCCAGGTTCAT--C0                    0               2
## AAAGTCCTCCGAGGCT--C0                    0               2
## AACAAAGCACCTGCGA--C0                    0               2
## AACCTGACATGCGGTC--C0                    0               3
## AACGTCAAGGTTCATC--C0                    0               7
##                      integrated_snn_res.0.1 integrated_snn_res.0.2
## AAAGGGCAGTCACGAG--C0                      0                      0
## AAAGTCCCAGGTTCAT--C0                      0                      0
## AAAGTCCTCCGAGGCT--C0                      0                      0
## AACAAAGCACCTGCGA--C0                      0                      0
## AACCTGACATGCGGTC--C0                      0                      0
## AACGTCAAGGTTCATC--C0                      1                      2
##                      integrated_snn_res.0.3 integrated_snn_res.0.4
## AAAGGGCAGTCACGAG--C0                      0                      0
## AAAGTCCCAGGTTCAT--C0                      0                      0
## AAAGTCCTCCGAGGCT--C0                      0                      0
## AACAAAGCACCTGCGA--C0                      0                      0
## AACCTGACATGCGGTC--C0                      0                      0
## AACGTCAAGGTTCATC--C0                      2                      2
##                      integrated_snn_res.0.5 integrated_snn_res.0.6
## AAAGGGCAGTCACGAG--C0                      0                      0
## AAAGTCCCAGGTTCAT--C0                      0                      0
## AAAGTCCTCCGAGGCT--C0                      0                      0
## AACAAAGCACCTGCGA--C0                      0                      0
## AACCTGACATGCGGTC--C0                      0                      0
## AACGTCAAGGTTCATC--C0                      2                      2
##                      integrated_snn_res.0.7 integrated_snn_res.0.8
## AAAGGGCAGTCACGAG--C0                      0                      0
## AAAGTCCCAGGTTCAT--C0                      0                      0
## AAAGTCCTCCGAGGCT--C0                      0                      0
## AACAAAGCACCTGCGA--C0                      0                      0
## AACCTGACATGCGGTC--C0                      0                      0
## AACGTCAAGGTTCATC--C0                      2                      2
##                      integrated_snn_res.0.9 integrated_snn_res.1
## AAAGGGCAGTCACGAG--C0                      0                    0
## AAAGTCCCAGGTTCAT--C0                      0                    0
## AAAGTCCTCCGAGGCT--C0                      0                    0
## AACAAAGCACCTGCGA--C0                      0                    0
## AACCTGACATGCGGTC--C0                      0                    0
## AACGTCAAGGTTCATC--C0                      2                    2
##                      integrated_snn_res.1.1 integrated_snn_res.1.2
## AAAGGGCAGTCACGAG--C0                      0                      3
## AAAGTCCCAGGTTCAT--C0                      0                      0
## AAAGTCCTCCGAGGCT--C0                      0                      0
## AACAAAGCACCTGCGA--C0                      0                      0
## AACCTGACATGCGGTC--C0                      0                      0
## AACGTCAAGGTTCATC--C0                      2                      2
##                      integrated_snn_res.1.3 integrated_snn_res.1.4
## AAAGGGCAGTCACGAG--C0                      2                      2
## AAAGTCCCAGGTTCAT--C0                      1                      1
## AAAGTCCTCCGAGGCT--C0                      1                      1
## AACAAAGCACCTGCGA--C0                      2                      2
## AACCTGACATGCGGTC--C0                      2                      2
## AACGTCAAGGTTCATC--C0                      3                      3
##                      integrated_snn_res.1.5 integrated_snn_res.1.6
## AAAGGGCAGTCACGAG--C0                      2                      3
## AAAGTCCCAGGTTCAT--C0                      0                      0
## AAAGTCCTCCGAGGCT--C0                      0                      4
## AACAAAGCACCTGCGA--C0                      2                      4
## AACCTGACATGCGGTC--C0                      2                      3
## AACGTCAAGGTTCATC--C0                      3                      2
##                      integrated_snn_res.1.7 integrated_snn_res.1.8
## AAAGGGCAGTCACGAG--C0                      3                      3
## AAAGTCCCAGGTTCAT--C0                      0                      2
## AAAGTCCTCCGAGGCT--C0                      4                      2
## AACAAAGCACCTGCGA--C0                      4                      2
## AACCTGACATGCGGTC--C0                      3                      3
## AACGTCAAGGTTCATC--C0                      1                      7
##                      integrated_snn_res.1.9 integrated_snn_res.2
## AAAGGGCAGTCACGAG--C0                      3                    6
## AAAGTCCCAGGTTCAT--C0                      2                    2
## AAAGTCCTCCGAGGCT--C0                      2                    2
## AACAAAGCACCTGCGA--C0                      2                    2
## AACCTGACATGCGGTC--C0                      8                    3
## AACGTCAAGGTTCATC--C0                      7                    7
##                      integrated_snn_res.2.1 integrated_snn_res.2.2
## AAAGGGCAGTCACGAG--C0                     26                     29
## AAAGTCCCAGGTTCAT--C0                     18                     22
## AAAGTCCTCCGAGGCT--C0                     18                     22
## AACAAAGCACCTGCGA--C0                     18                     22
## AACCTGACATGCGGTC--C0                     26                     29
## AACGTCAAGGTTCATC--C0                     36                     40
##                      integrated_snn_res.2.3 integrated_snn_res.2.4
## AAAGGGCAGTCACGAG--C0                     31                     32
## AAAGTCCCAGGTTCAT--C0                     22                     20
## AAAGTCCTCCGAGGCT--C0                     22                     20
## AACAAAGCACCTGCGA--C0                     22                     20
## AACCTGACATGCGGTC--C0                     31                     32
## AACGTCAAGGTTCATC--C0                     40                     41
##                      integrated_snn_res.2.5 integrated_snn_res.2.6
## AAAGGGCAGTCACGAG--C0                     30                     31
## AAAGTCCCAGGTTCAT--C0                     23                     20
## AAAGTCCTCCGAGGCT--C0                     23                     20
## AACAAAGCACCTGCGA--C0                     23                     20
## AACCTGACATGCGGTC--C0                     30                     31
## AACGTCAAGGTTCATC--C0                     41                     42
##                      integrated_snn_res.2.7 integrated_snn_res.2.8
## AAAGGGCAGTCACGAG--C0                     33                     33
## AAAGTCCCAGGTTCAT--C0                     23                     22
## AAAGTCCTCCGAGGCT--C0                     23                     22
## AACAAAGCACCTGCGA--C0                     23                     22
## AACCTGACATGCGGTC--C0                     33                     33
## AACGTCAAGGTTCATC--C0                     44                     44
##                      integrated_snn_res.2.9 integrated_snn_res.3
## AAAGGGCAGTCACGAG--C0                     32                   37
## AAAGTCCCAGGTTCAT--C0                     22                   22
## AAAGTCCTCCGAGGCT--C0                     22                   22
## AACAAAGCACCTGCGA--C0                     22                   22
## AACCTGACATGCGGTC--C0                     32                   37
## AACGTCAAGGTTCATC--C0                     43                   46
##                      integrated_snn_res.0.21 integrated_snn_res.0.22
## AAAGGGCAGTCACGAG--C0                       6                       6
## AAAGTCCCAGGTTCAT--C0                       6                       6
## AAAGTCCTCCGAGGCT--C0                       6                       6
## AACAAAGCACCTGCGA--C0                       6                       6
## AACCTGACATGCGGTC--C0                       6                       6
## AACGTCAAGGTTCATC--C0                       6                       6
##                      integrated_snn_res.0.23 integrated_snn_res.0.24
## AAAGGGCAGTCACGAG--C0                       6                       7
## AAAGTCCCAGGTTCAT--C0                       6                       7
## AAAGTCCTCCGAGGCT--C0                       6                       7
## AACAAAGCACCTGCGA--C0                       6                       7
## AACCTGACATGCGGTC--C0                       6                       7
## AACGTCAAGGTTCATC--C0                       6                       7
##                      integrated_snn_res.0.25 integrated_snn_res.0.26
## AAAGGGCAGTCACGAG--C0                       6                       6
## AAAGTCCCAGGTTCAT--C0                       6                       6
## AAAGTCCTCCGAGGCT--C0                       6                       6
## AACAAAGCACCTGCGA--C0                       6                       6
## AACCTGACATGCGGTC--C0                       6                       6
## AACGTCAAGGTTCATC--C0                       6                       6
##                      integrated_snn_res.0.27 integrated_snn_res.0.28
## AAAGGGCAGTCACGAG--C0                       6                       6
## AAAGTCCCAGGTTCAT--C0                       6                       6
## AAAGTCCTCCGAGGCT--C0                       6                       6
## AACAAAGCACCTGCGA--C0                       6                       6
## AACCTGACATGCGGTC--C0                       6                       6
## AACGTCAAGGTTCATC--C0                       6                       6
##                      integrated_snn_res.0.29 annotation.1 celltype.stim
## AAAGGGCAGTCACGAG--C0                       6        Macro      Macro_C0
## AAAGTCCCAGGTTCAT--C0                       6        Macro      Macro_C0
## AAAGTCCTCCGAGGCT--C0                       6        Macro      Macro_C0
## AACAAAGCACCTGCGA--C0                       6        Macro      Macro_C0
## AACCTGACATGCGGTC--C0                       6        Macro      Macro_C0
## AACGTCAAGGTTCATC--C0                       6        Macro      Macro_C0
##                      celltype res.0.5.stim annotated.1.sample.id
## AAAGGGCAGTCACGAG--C0    Macro         7_C0              Macro_C0
## AAAGTCCCAGGTTCAT--C0    Macro         7_C0              Macro_C0
## AAAGTCCTCCGAGGCT--C0    Macro         7_C0              Macro_C0
## AACAAAGCACCTGCGA--C0    Macro         7_C0              Macro_C0
## AACCTGACATGCGGTC--C0    Macro         7_C0              Macro_C0
## AACGTCAAGGTTCATC--C0    Macro         7_C0              Macro_C0
##                      macrophages_subclusters_res.0
## AAAGGGCAGTCACGAG--C0                             0
## AAAGTCCCAGGTTCAT--C0                             0
## AAAGTCCTCCGAGGCT--C0                             0
## AACAAAGCACCTGCGA--C0                             0
## AACCTGACATGCGGTC--C0                             0
## AACGTCAAGGTTCATC--C0                             0
##                      macrophages_subclusters_res.0.1
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               1
##                      macrophages_subclusters_res.0.2
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.0.3
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.0.4
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.0.5
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.0.6
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.0.7
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.0.8
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.0.9
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.1
## AAAGGGCAGTCACGAG--C0                             0
## AAAGTCCCAGGTTCAT--C0                             0
## AAAGTCCTCCGAGGCT--C0                             0
## AACAAAGCACCTGCGA--C0                             0
## AACCTGACATGCGGTC--C0                             0
## AACGTCAAGGTTCATC--C0                             2
##                      macrophages_subclusters_res.1.1
## AAAGGGCAGTCACGAG--C0                               0
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.1.2
## AAAGGGCAGTCACGAG--C0                               3
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               0
## AACCTGACATGCGGTC--C0                               0
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.1.3
## AAAGGGCAGTCACGAG--C0                               2
## AAAGTCCCAGGTTCAT--C0                               1
## AAAGTCCTCCGAGGCT--C0                               1
## AACAAAGCACCTGCGA--C0                               2
## AACCTGACATGCGGTC--C0                               2
## AACGTCAAGGTTCATC--C0                               3
##                      macrophages_subclusters_res.1.4
## AAAGGGCAGTCACGAG--C0                               2
## AAAGTCCCAGGTTCAT--C0                               1
## AAAGTCCTCCGAGGCT--C0                               1
## AACAAAGCACCTGCGA--C0                               2
## AACCTGACATGCGGTC--C0                               2
## AACGTCAAGGTTCATC--C0                               3
##                      macrophages_subclusters_res.1.5
## AAAGGGCAGTCACGAG--C0                               2
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               0
## AACAAAGCACCTGCGA--C0                               2
## AACCTGACATGCGGTC--C0                               2
## AACGTCAAGGTTCATC--C0                               3
##                      macrophages_subclusters_res.1.6
## AAAGGGCAGTCACGAG--C0                               3
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               4
## AACAAAGCACCTGCGA--C0                               4
## AACCTGACATGCGGTC--C0                               3
## AACGTCAAGGTTCATC--C0                               2
##                      macrophages_subclusters_res.1.7
## AAAGGGCAGTCACGAG--C0                               3
## AAAGTCCCAGGTTCAT--C0                               0
## AAAGTCCTCCGAGGCT--C0                               4
## AACAAAGCACCTGCGA--C0                               4
## AACCTGACATGCGGTC--C0                               3
## AACGTCAAGGTTCATC--C0                               1
##                      macrophages_subclusters_res.1.8
## AAAGGGCAGTCACGAG--C0                               3
## AAAGTCCCAGGTTCAT--C0                               2
## AAAGTCCTCCGAGGCT--C0                               2
## AACAAAGCACCTGCGA--C0                               2
## AACCTGACATGCGGTC--C0                               3
## AACGTCAAGGTTCATC--C0                               7
##                      macrophages_subclusters_res.1.9
## AAAGGGCAGTCACGAG--C0                               3
## AAAGTCCCAGGTTCAT--C0                               2
## AAAGTCCTCCGAGGCT--C0                               2
## AACAAAGCACCTGCGA--C0                               2
## AACCTGACATGCGGTC--C0                               8
## AACGTCAAGGTTCATC--C0                               7
##                      macrophages_subclusters_res.2
## AAAGGGCAGTCACGAG--C0                             6
## AAAGTCCCAGGTTCAT--C0                             2
## AAAGTCCTCCGAGGCT--C0                             2
## AACAAAGCACCTGCGA--C0                             2
## AACCTGACATGCGGTC--C0                             3
## AACGTCAAGGTTCATC--C0                             7
# install.packages('clustree')
library(clustree)
c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.")
c

pdf(paste0(outdir, "/2_clustree_macrophages.pdf"))
c
dev.off()
## png 
##   2

Using the stability index to assess clusters The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions and is automatically calculated when a clustering tree is built.

c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.", 
    node_colour = "sc3_stability")
pdf(paste0(outdir, "/3_clustree_macrophages_sc3_stability.pdf"))
c
dev.off()
## png 
##   2

SOCS3, IL-1R2, IL1rn

Socs3, Il1r2, Il1rn

c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.", 
    node_colour = "Socs3", node_colour_aggr = "mean")
c

pdf(paste0(outdir, "/4_clustree_macrophages_Socs3.pdf"))
c
dev.off()
## png 
##   2
c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.", 
    node_colour = "Il1r2", node_colour_aggr = "mean")
c

pdf(paste0(outdir, "/5_clustree_macrophages_Il1r2.pdf"))
c
dev.off()
## png 
##   2
c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.", 
    node_colour = "Il1rn", node_colour_aggr = "mean")
c

pdf(paste0(outdir, "/6_clustree_macrophages_Il1rn.pdf"))
c
dev.off()
## png 
##   2

RES 0.3

# install.packages('Polychrome')
library(Polychrome)
set.seed(5658807)  # for reproducibility
levels <- levels(as.factor(macrophages.integrated$macrophages_subclusters_res.0.3))
rainbow2.6c <- c("#03539C", "#12999E", "#B7CE05", "#FAEB09", 
    "#FDA908", "#E82564")
colors.macrophages.subclusters <- createPalette(length(levels), 
    rainbow2.6c, M = 1000)
names(colors.macrophages.subclusters) <- levels
colors.macrophages.subclusters
##         0         1         2         3      <NA>      <NA>      <NA>      <NA> 
## "#1662B5" "#0D9A9D" "#C2DB1C" "#FBE22E" "#FCA700" "#F33877" "#FA35F5" "#FFB1DB" 
##      <NA>      <NA> 
## "#40F897" "#63691C"
slices <- rep(1, length(colors.macrophages.subclusters))
pie(slices, col = colors.macrophages.subclusters, labels = names(colors.macrophages.subclusters))

Idents(macrophages.integrated) <- "macrophages_subclusters_res.0.3"
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 8, cols = colors.macrophages.subclusters)

pdf(paste0(outdir, "/7_DimPlot_umap_macrophages_subclusters_res.0.3.pdf"))
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 0.5, 
    label = T, label.size = 6, cols = colors.macrophages.subclusters)
dev.off()
## png 
##   2
colors.samples <- c("#12999E", "#FDA908")
names(colors.samples) <- c("C0", "C24")
macrophages.integrated$sample.id <- factor(macrophages.integrated$sample.id, 
    levels = names(colors.samples))
p <- DimPlot(macrophages.integrated, reduction = "umap", group.by = "sample.id", 
    cols = colors.samples, label = F)
pdf(paste0(outdir, "/8_UMAP_macrophages_samples.pdf"), width = 10, 
    height = 9)
p
dev.off()
## png 
##   2
pdf(paste0(outdir, "/9_UMAP_macrophages_samples_noLegend.pdf"), 
    width = 7, height = 7)
p + NoLegend()
dev.off()
## png 
##   2

Rename subclusters

new.cluster.ids <- paste0("Macro", 1:length(levels(as.factor(macrophages.integrated$macrophages_subclusters_res.0.3))))
names(new.cluster.ids) <- levels(macrophages.integrated)
macrophages.integrated <- RenameIdents(macrophages.integrated, 
    new.cluster.ids)
macrophages.integrated[["renamed.macrophages.subclusters"]] <- macrophages.integrated@active.ident
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
names(colors.macrophages.subclusters) <- levels(macrophages.integrated$renamed.macrophages.subclusters)
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 8, cols = colors.macrophages.subclusters)

pdf(paste0(outdir, "/10_DimPlot_umap_renamed_macrophages_subclusters.pdf"))
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 0.5, 
    label = T, label.size = 6, cols = colors.macrophages.subclusters)
dev.off()
## png 
##   2
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 1, 
    label = F, label.size = 8, cols = colors.macrophages.subclusters) + 
    NoLegend()

pdf(paste0(outdir, "/11_DimPlot_umap_renamed_macrophages_subclusters_noLabels_noLegend.pdf"))
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 0.5, 
    label = F, label.size = 6, cols = colors.macrophages.subclusters) + 
    NoLegend()
dev.off()
## png 
##   2
saveRDS(macrophages.integrated, paste0(outdir, "/12.macrophages.integrated.rds"))

Conserved markers

t <- table(macrophages.integrated$renamed.macrophages.subclusters, 
    macrophages.integrated$sample.id)
t
##         
##           C0 C24
##   Macro1 137 161
##   Macro2  72  98
##   Macro3  39  50
##   Macro4  19  12
openxlsx::write.xlsx(as.data.frame.matrix(t), file = paste0(outdir, 
    "/13_Sample_macrophages_subcluster_composition.xlsx"), rowNames = T, 
    colNames = T)
DefaultAssay(macrophages.integrated) <- "RNA"
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
subclusters <- levels(as.factor(macrophages.integrated$renamed.macrophages.subclusters))
conserved.markers <- data.frame(matrix(ncol = 14))
for (c in subclusters) {
    print(c)
    markers.c <- FindConservedMarkers(macrophages.integrated, 
        ident.1 = c, grouping.var = "sample.id", verbose = T, 
        logfc.threshold = -Inf, min.pct = -Inf, min.cells.feature = -Inf, 
        min.cells.group = -Inf)
    markers.c <- cbind(data.frame(cluster = rep(c, dim(markers.c)[1]), 
        gene = rownames(markers.c)), markers.c)
    write.table(markers.c, file = paste0(outdir, "/14_markers_", 
        c, "_macrophages_subclusters.txt"))
    colnames(conserved.markers) <- colnames(markers.c)
    conserved.markers <- rbind(conserved.markers, markers.c)
}
## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"
## [1] "Macro4"
# for (c in 1:4) { subcluster <- paste0('AM', c) markers.c <-
# read.table(file = paste0(outdir, '/31_markers_',
# subcluster, '_AM_subclusters.txt'))
# colnames(conserved.markers) <- colnames(markers.c)
# conserved.markers <- rbind(conserved.markers, markers.c) }
# head(conserved.markers)
# markers.AM5 <- read.table(file = paste0(outdir,
# '/31_markers_', 'AM5', '_AM_subclusters.txt'))
# colnames(markers.AM5) colnames(conserved.markers)
# for ( i in 1:10 ) { markers.AM5[, 34+i] <- rep(NA,
# dim(markers.AM5)[1]) } markers.AM5 <- markers.AM5[, c(1:12,
# 35:44, 13:34)] colnames(markers.AM5) <-
# colnames(conserved.markers) head(markers.AM5)
# conserved.markers <- rbind(conserved.markers, markers.AM5)
conserved.markers <- conserved.markers[-1, ]
conserved.markers <- conserved.markers[, c(1:2, 13:14, 3:12)]
head(conserved.markers)
##       cluster  gene     max_pval minimump_p_val     C0_p_val C0_avg_log2FC
## Sdc4   Macro1  Sdc4 6.664112e-02   1.165391e-30 6.664112e-02    0.14916480
## Cd14   Macro1  Cd14 6.988969e-19   5.627031e-30 6.988969e-19    1.71602690
## Cebpb  Macro1 Cebpb 1.710274e-03   2.737600e-24 1.710274e-03    0.81625576
## C1qc   Macro1  C1qc 4.471371e-17   4.853862e-23 4.471371e-17    1.15050250
## Socs3  Macro1 Socs3 6.343829e-01   3.133599e-22 6.343829e-01    0.09794724
## Pid1   Macro1  Pid1 2.189482e-02   3.136037e-22 2.189482e-02    0.23031958
##       C0_pct.1 C0_pct.2 C0_p_val_adj    C24_p_val C24_avg_log2FC C24_pct.1
## Sdc4     0.328    0.223 1.000000e+00 5.826955e-31       1.547758     0.994
## Cd14     0.803    0.277 1.565459e-14 2.813516e-30       2.000407     0.981
## Cebpb    0.219    0.092 1.000000e+00 1.368800e-24       1.445991     0.988
## C1qc     0.964    0.546 1.001542e-12 2.426931e-23       1.415284     0.957
## Socs3    0.117    0.100 1.000000e+00 1.566800e-22       1.329294     0.938
## Pid1     0.453    0.331 1.000000e+00 1.568018e-22       1.396260     0.776
##       C24_pct.2 C24_p_val_adj
## Sdc4      0.656  1.305180e-26
## Cd14      0.475  6.301994e-26
## Cebpb     0.537  3.065975e-20
## C1qc      0.369  5.436082e-19
## Socs3     0.469  3.509474e-18
## Pid1      0.225  3.512205e-18

Only top markers that are positive markers

conserved.markers$marker.type <- apply(conserved.markers, 1, function(x) {
  y <- as.numeric(x)
  if ( (if (is.na(y[6])) {TRUE} else {y[6]>0})
       & (if (is.na(y[11])) {TRUE} else {y[11]>0})
       # & (if (is.na(y[16])) {TRUE} else {y[16]>0})
       # & (if (is.na(y[21])) {TRUE} else {y[21]>0})
       # & (if (is.na(y[26])) {TRUE} else {y[26]>0})
       # & (if (is.na(y[31])) {TRUE} else {y[31]>0})
       # & (if (is.na(y[36])) {TRUE} else {y[36]>0})
       # & (if (is.na(y[41])) {TRUE} else {y[41]>0})
       )
    {"POS"}
  else if ( ( if (is.na(y[6])) {TRUE} else {y[6]<0})
       & (if (is.na(y[11])) {TRUE} else {y[11]<0})
       # & (if (is.na(y[16])) {TRUE} else {y[16]<0})
       # & (if (is.na(y[21])) {TRUE} else {y[21]<0})
       # & (if (is.na(y[26])) {TRUE} else {y[26]<0})
       # & (if (is.na(y[31])) {TRUE} else {y[31]<0})
       # & (if (is.na(y[36])) {TRUE} else {y[36]<0})
       # & (if (is.na(y[41])) {TRUE} else {y[41]<0})
       ) 
    {"NEG"}
  else {"UND"}
  })
conserved.markers <- conserved.markers[, c(1:4, 15, 5:14)]
openxlsx::write.xlsx(conserved.markers, paste0(outdir, "/15_conserved_markers_macrophages.xlsx"), 
    colNames = T)
head(conserved.markers)
##       cluster  gene     max_pval minimump_p_val marker.type     C0_p_val
## Sdc4   Macro1  Sdc4 6.664112e-02   1.165391e-30         POS 6.664112e-02
## Cd14   Macro1  Cd14 6.988969e-19   5.627031e-30         POS 6.988969e-19
## Cebpb  Macro1 Cebpb 1.710274e-03   2.737600e-24         POS 1.710274e-03
## C1qc   Macro1  C1qc 4.471371e-17   4.853862e-23         POS 4.471371e-17
## Socs3  Macro1 Socs3 6.343829e-01   3.133599e-22         POS 6.343829e-01
## Pid1   Macro1  Pid1 2.189482e-02   3.136037e-22         POS 2.189482e-02
##       C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj    C24_p_val C24_avg_log2FC
## Sdc4     0.14916480    0.328    0.223 1.000000e+00 5.826955e-31       1.547758
## Cd14     1.71602690    0.803    0.277 1.565459e-14 2.813516e-30       2.000407
## Cebpb    0.81625576    0.219    0.092 1.000000e+00 1.368800e-24       1.445991
## C1qc     1.15050250    0.964    0.546 1.001542e-12 2.426931e-23       1.415284
## Socs3    0.09794724    0.117    0.100 1.000000e+00 1.566800e-22       1.329294
## Pid1     0.23031958    0.453    0.331 1.000000e+00 1.568018e-22       1.396260
##       C24_pct.1 C24_pct.2 C24_p_val_adj
## Sdc4      0.994     0.656  1.305180e-26
## Cd14      0.981     0.475  6.301994e-26
## Cebpb     0.988     0.537  3.065975e-20
## C1qc      0.957     0.369  5.436082e-19
## Socs3     0.938     0.469  3.509474e-18
## Pid1      0.776     0.225  3.512205e-18
saveRDS(macrophages.integrated, paste0(outdir, "/16.macrophages.integrated.rds"))
conserved.markers.pos <- conserved.markers[conserved.markers$marker.type == 
    "POS", ]
head(conserved.markers.pos)
##       cluster  gene     max_pval minimump_p_val marker.type     C0_p_val
## Sdc4   Macro1  Sdc4 6.664112e-02   1.165391e-30         POS 6.664112e-02
## Cd14   Macro1  Cd14 6.988969e-19   5.627031e-30         POS 6.988969e-19
## Cebpb  Macro1 Cebpb 1.710274e-03   2.737600e-24         POS 1.710274e-03
## C1qc   Macro1  C1qc 4.471371e-17   4.853862e-23         POS 4.471371e-17
## Socs3  Macro1 Socs3 6.343829e-01   3.133599e-22         POS 6.343829e-01
## Pid1   Macro1  Pid1 2.189482e-02   3.136037e-22         POS 2.189482e-02
##       C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj    C24_p_val C24_avg_log2FC
## Sdc4     0.14916480    0.328    0.223 1.000000e+00 5.826955e-31       1.547758
## Cd14     1.71602690    0.803    0.277 1.565459e-14 2.813516e-30       2.000407
## Cebpb    0.81625576    0.219    0.092 1.000000e+00 1.368800e-24       1.445991
## C1qc     1.15050250    0.964    0.546 1.001542e-12 2.426931e-23       1.415284
## Socs3    0.09794724    0.117    0.100 1.000000e+00 1.566800e-22       1.329294
## Pid1     0.23031958    0.453    0.331 1.000000e+00 1.568018e-22       1.396260
##       C24_pct.1 C24_pct.2 C24_p_val_adj
## Sdc4      0.994     0.656  1.305180e-26
## Cd14      0.981     0.475  6.301994e-26
## Cebpb     0.988     0.537  3.065975e-20
## C1qc      0.957     0.369  5.436082e-19
## Socs3     0.938     0.469  3.509474e-18
## Pid1      0.776     0.225  3.512205e-18
table(macrophages.integrated$sample.id)
## 
##  C0 C24 
## 267 321
colnames(conserved.markers.pos)
##  [1] "cluster"        "gene"           "max_pval"       "minimump_p_val"
##  [5] "marker.type"    "C0_p_val"       "C0_avg_log2FC"  "C0_pct.1"      
##  [9] "C0_pct.2"       "C0_p_val_adj"   "C24_p_val"      "C24_avg_log2FC"
## [13] "C24_pct.1"      "C24_pct.2"      "C24_p_val_adj"
library(dplyr)
top1.pos <- conserved.markers.pos %>% group_by(cluster) %>% top_n(n = 1, 
    wt = -minimump_p_val)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = -max_pval)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = C24_avg_log2FC)
head(top1.pos)
## # A tibble: 4 x 15
## # Groups:   cluster [4]
##   cluster gene  max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
##   <chr>   <chr>    <dbl>          <dbl> <chr>          <dbl>         <dbl>
## 1 Macro1  Sdc4  6.66e- 2       1.17e-30 POS         6.66e- 2         0.149
## 2 Macro2  Kap   3.17e-12       1.10e-21 POS         3.17e-12         0.748
## 3 Macro3  Pfkp  2.36e- 2       4.23e-39 POS         2.36e- 2         0.176
## 4 Macro4  Rets… 5.04e- 5       2.83e-21 POS         1.42e-21         0.998
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## #   C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## #   C24_p_val_adj <dbl>
DefaultAssay(macrophages.integrated) <- "RNA"
f <- FeaturePlot(macrophages.integrated, features = top1.pos$gene, 
    min.cutoff = "q9", order = T, label = T, label.size = 8, 
    repel = T)
f

pdf(paste0(outdir, "/17_FeaturePlot_top_pos_markers_macrophages.pdf"), 
    width = 14, height = 14)
f
dev.off()
## png 
##   2

DotPlot of top markers

DefaultAssay(macrophages.integrated) <- "RNA"
d <- DotPlot(object = macrophages.integrated, features = top1.pos$gene, 
    cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/18_DotPlot_top_cell_types_pos_markers_clusters.macrophages.pdf"), 
    width = 7, height = 7)
d
dev.off()
## png 
##   2

DotPlot of top 10 markers

library(dplyr)
top10.pos <- conserved.markers.pos %>% group_by(cluster) %>% 
    top_n(n = 10, wt = -minimump_p_val)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10, 
    wt = -max_pval)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10, 
    wt = C24_avg_log2FC)
head(top10.pos)
## # A tibble: 6 x 15
## # Groups:   cluster [1]
##   cluster gene  max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
##   <chr>   <chr>    <dbl>          <dbl> <chr>          <dbl>         <dbl>
## 1 Macro1  Sdc4  6.66e- 2       1.17e-30 POS         6.66e- 2        0.149 
## 2 Macro1  Cd14  6.99e-19       5.63e-30 POS         6.99e-19        1.72  
## 3 Macro1  Cebpb 1.71e- 3       2.74e-24 POS         1.71e- 3        0.816 
## 4 Macro1  C1qc  4.47e-17       4.85e-23 POS         4.47e-17        1.15  
## 5 Macro1  Socs3 6.34e- 1       3.13e-22 POS         6.34e- 1        0.0979
## 6 Macro1  Pid1  2.19e- 2       3.14e-22 POS         2.19e- 2        0.230 
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## #   C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## #   C24_p_val_adj <dbl>
DefaultAssay(macrophages.integrated) <- "RNA"
d <- DotPlot(object = macrophages.integrated, features = top10.pos$gene, 
    cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/19_DotPlot_top10_cell_types_pos_markers_clusters.macrophages.pdf"), 
    width = 14, height = 7)
d
dev.off()
## png 
##   2

Segregation of clusters by various sources of uninteresting variation

Explore additional metrics, such as the number of UMIs and genes per cell, and mitochondrial gene expression by UMAP.

# Determine metrics to plot present in
# tcells.integrated@meta.data
metrics <- c("nCount_RNA", "nFeature_RNA", "percent.mito")

f <- FeaturePlot(macrophages.integrated, reduction = "umap", 
    features = metrics, pt.size = 0.4, order = TRUE, min.cutoff = "q10", 
    label = TRUE)
f

pdf(paste0(outdir, "/20_FeaturePlot_macrophages.integrated_umap_metrics.pdf"), 
    width = 14, height = 14)
f
dev.off()
## png 
##   2

Heatmap

macrophages.integrated <- ScaleData(macrophages.integrated, verbose = T, 
    features = rownames(macrophages.integrated))
d <- DoHeatmap(macrophages.integrated, features = top10.pos$gene, 
    group.by = "renamed.macrophages.subclusters", group.colors = colors.macrophages.subclusters) + 
    NoLegend()
d

pdf(paste0(outdir, "/21_DoHeatmap_top10_pos_macrophages.pdf"), 
    width = 15, height = 13)
d
dev.off()
## png 
##   2

DotPlot with known/useful markers

immune.cells <- c("Ptprc", "Spi1", "Cd3g")
b.cells <- c("Cd19", "Cd79a", "Cd79b")
t.cells <- "Cd3e"
nk.cells <- c("Ncr1")
myeloid.cells <- "Cd68"
neutrophils <- c("Csf3r", "S100a8", "S100a9", "Retnlg")
monocytes.macrophages.dendritic.cells <- "Csf1r"
dc <- c("Siglech", "Cd209a")
dc1 <- "Xcr1"
dc2 <- c("Ccl18", "Nr4a3", "Tnfsf12")
dc3 <- "Fscn1"
pdc <- "Tcf4"
monocytes <- c("Ly6c1", "Ly6c2", "Ccr2", "Cx3cr1", "Fcgr4")
macrophages <- c("Mrc1", "Lyz2", "Adgre1", "Axl", "Cxcl2", "Trem2", 
    "C1qc", "Fcgr1", "C5ar1")
am.markers <- c("Pparg", "Itgax", "Ear1", "Ear2", "Ear3", "Car4", 
    "Siglec5", "SiglecF")
im.markers <- c("Itgam", "Cx3cr1", "Ccl12", "Bhlme40", "Trem2", 
    "Cxcl12", "Csf1")
peribronchial.macrophages <- c("Ccr2", "Bhlhe40", "Bhlhe41")
perivascular.macrophages <- c("Lyve1", "F13a1")
proliferating.cells <- "Mki67"
mesenchymal.cells <- c("Col1a1", "Col1a2", "Pdgfra")
epithelial.cells <- c("Epcam", "Cdh1")
mast.cells <- "Mcpt8"
endothelial.cells <- c("Pecam1", "Emcn")
markers <- unique(c(immune.cells, b.cells, t.cells, nk.cells, 
    myeloid.cells, neutrophils, monocytes.macrophages.dendritic.cells, 
    dc, dc1, dc2, dc3, pdc, monocytes, macrophages, am.markers, 
    im.markers, peribronchial.macrophages, perivascular.macrophages, 
    proliferating.cells, mesenchymal.cells, epithelial.cells, 
    mast.cells, endothelial.cells))

Dot plot

DefaultAssay(macrophages.integrated) <- "RNA"
d <- DotPlot(object = macrophages.integrated, features = markers, 
    cols = c("#03539C", "#E82564"), dot.scale = 15, group.by = "renamed.macrophages.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/22_DotPlot_known_markers_macrophages.pdf"), 
    width = 18, height = 4)
d
dev.off()
## png 
##   2

Expression of genes of interest: Socs3, Il1r2, Il1rn

goi <- c("Socs3", "Il1r2", "Il1rn")
DefaultAssay(macrophages.integrated) <- "RNA"
f <- FeaturePlot(macrophages.integrated, features = goi, min.cutoff = "q9", 
    order = T, label = T, label.size = 8, repel = T)
f

pdf(paste0(outdir, "/23_FeaturePlot_goi_macrophages.pdf"), width = 14, 
    height = 14)
f
dev.off()
## png 
##   2
DefaultAssay(macrophages.integrated) <- "RNA"
d <- DotPlot(object = macrophages.integrated, features = goi, 
    cols = c("#03539C", "#E82564"), dot.scale = 15, group.by = "renamed.macrophages.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/24_DotPlot_goi_macrophages.pdf"), width = 6, 
    height = 4)
d
dev.off()
## png 
##   2

DIFFERENTIAL EXPRESSION ANALYSES

differential expression studies for the following groups: 1. C24 vs C0 within each cell type (from manual annotation annotation.1)

1. C24 vs C0 within each cell type (from manual annotation annotation.1)

macrophages.integrated$renamed.macrophages.subclusters.sample.id <- paste(macrophages.integrated$renamed.macrophages.subclusters, 
    macrophages.integrated$sample.id, sep = "_")
levels(factor(macrophages.integrated$renamed.macrophages.subclusters.sample.id))
## [1] "Macro1_C0"  "Macro1_C24" "Macro2_C0"  "Macro2_C24" "Macro3_C0" 
## [6] "Macro3_C24" "Macro4_C0"  "Macro4_C24"
DefaultAssay(macrophages.integrated) <- "RNA"
output <- paste0(outdir, "/25_DE_macrophages_subclusters_C24_vs_C0_")
levels(factor(macrophages.integrated$renamed.macrophages.subclusters.sample.id))
## [1] "Macro1_C0"  "Macro1_C24" "Macro2_C0"  "Macro2_C24" "Macro3_C0" 
## [6] "Macro3_C24" "Macro4_C0"  "Macro4_C24"
for (cell.type in levels(factor(macrophages.integrated$renamed.macrophages.subclusters))) {
    try({
        print(cell.type)
        ident1 <- paste0(cell.type, "_C24")
        ident2 <- paste0(cell.type, "_C0")
        Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
        condition.diffgenes <- FindMarkers(macrophages.integrated, 
            ident.1 = ident1, ident.2 = ident2, logfc.threshold = -Inf, 
            test.use = "wilcox", min.pct = -Inf, verbose = T, 
            min.cells.feature = -Inf, min.cells.group = -Inf)
        
        condition.diffgenes.adjpval.0.05 <- condition.diffgenes[condition.diffgenes$p_val_adj < 
            0.05, ]
        condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05[condition.diffgenes.adjpval.0.05$avg_log2FC > 
            0, ]
        condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05.upregulated.sorted[order(condition.diffgenes.adjpval.0.05.upregulated.sorted$avg_log2FC, 
            decreasing = T), ]
        condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05[condition.diffgenes.adjpval.0.05$avg_log2FC < 
            0, ]
        condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05.downregulated.sorted[order(condition.diffgenes.adjpval.0.05.downregulated.sorted$avg_log2FC, 
            decreasing = F), ]
        list_of_datasets <- list(sig.upregulated.genes.ordered.by.avg.logFC = condition.diffgenes.adjpval.0.05.upregulated.sorted, 
            sig.downregulated.genes.ordered.by.avg.logFC = condition.diffgenes.adjpval.0.05.downregulated.sorted, 
            all.genes.ordered.by.adj.pval = condition.diffgenes)
        openxlsx::write.xlsx(list_of_datasets, file = paste0(output, 
            cell.type, ".xlsx"), row.names = T)
    })
}
## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"
## [1] "Macro4"
for (cell.type in levels(factor(macrophages.integrated$renamed.macrophages.subclusters))) {
    try({
        print(cell.type)
        filename <- paste0(output, cell.type, ".xlsx")
        sheets <- openxlsx::getSheetNames(filename)
        SheetList <- lapply(sheets, openxlsx::read.xlsx, xlsxFile = filename)
        names(SheetList) <- sheets
        
        condition.diffgenes.adjpval.0.05.upregulated.sorted <- SheetList[[1]]
        rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted) <- condition.diffgenes.adjpval.0.05.upregulated.sorted[, 
            1]
        condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05.upregulated.sorted[, 
            -1]
        
        condition.diffgenes.adjpval.0.05.downregulated.sorted <- SheetList[[2]]
        rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted) <- condition.diffgenes.adjpval.0.05.downregulated.sorted[, 
            1]
        condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05.downregulated.sorted[, 
            -1]
        
        
        Idents(annotated) <- "renamed.macrophages.subclusters"
        
        try({
            f <- FeaturePlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1], 
                split.by = "sample.id", cols = c("grey", "#E82564"), 
                order = T, min.cutoff = "q10", max.cutoff = "q90")
            pdf(paste0(output, cell.type, "_FeaturePlot_top_upregulated_split_by_sample.pdf"), 
                width = 14, height = 7)
            print(f)
            dev.off()
        })
        
        try({
            f <- FeaturePlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1], 
                split.by = "sample.id", cols = c("grey", "#E82564"), 
                order = T, min.cutoff = "q10", max.cutoff = "q90")
            pdf(paste0(output, cell.type, "_FeaturePlot_top_downregulated_split_by_sample.pdf"), 
                width = 14, height = 7)
            print(f)
            dev.off()
        })
        
        try({
            plot <- VlnPlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1], 
                group.by = "renamed.macrophages.subclusters", 
                split.by = "sample.id", pt.size = 0, combine = FALSE, 
                cols = colors.samples)
            pdf(paste0(output, cell.type, "_VlnPlot_top_upregulated_group_by_cluster_split_by_sample.pdf"), 
                width = 11, height = 8.5)
            print(plot)
            dev.off()
        })
        
        try({
            plot <- VlnPlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1], 
                group.by = "renamed.macrophages.subclusters", 
                split.by = "sample.id", pt.size = 0, combine = FALSE, 
                cols = colors.samples)
            pdf(paste0(output, cell.type, "_VlnPlot_top_downregulated_group_by_cluster_split_by_sample.pdf"), 
                width = 11, height = 8.5)
            print(plot)
            dev.off()
        })
        
        ## DotPlot with 50 most upregulated genes
        Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
        try({
            d <- DotPlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1:50], 
                cols = c("#03539C", "#E82564"), dot.scale = 8, 
                group.by = "renamed.macrophages.subclusters.sample.id") + 
                RotatedAxis()
            pdf(paste0(output, cell.type, "_DotPlot_top50_upregulated_groupby_annotation.1.sample.id.pdf"), 
                width = (4 + min(50, length(rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted))) * 
                  15/50), height = 17)
            print(d)
            dev.off()
        })
        
        ## DotPlot with 50 most downregulated genes
        Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
        try({
            d <- DotPlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1:50], 
                cols = c("#03539C", "#E82564"), dot.scale = 8, 
                group.by = "renamed.macrophages.subclusters.sample.id") + 
                RotatedAxis()
            pdf(paste0(output, cell.type, "_DotPlot_top50_downregulated_groupby_annotation.1.sample.id.pdf"), 
                width = (4 + min(50, length(rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted))) * 
                  15/50), height = 17)
            print(d)
            dev.off()
        })
    })
}
## [1] "Macro1"
## [[1]]
## 
## [[1]]
## [1] "Macro2"
## [[1]]
## 
## [[1]]
## [1] "Macro3"
## [[1]]
## 
## [[1]]
## [1] "Macro4"
## Error : None of the requested features were found: NA in slot data
## [[1]]
## 
## Error in FetchData(object = object, vars = features, slot = slot) : 
##   None of the requested variables were found: NA
## Error in h(simpleError(msg, call)) : 
##   error in evaluating the argument 'x' in selecting a method for function 'mean': non-numeric argument to mathematical function

Plots for genes of interest

for (gene in goi) {
    try({
        Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
        try({
            f <- FeaturePlot(macrophages.integrated, features = gene, 
                split.by = "sample.id", cols = c("grey", "#E82564"), 
                order = T, min.cutoff = "q10", max.cutoff = "q90", 
                label = T)
            pdf(paste0(outdir, "/26_", gene, "_FeaturePlot_split_by_sample.pdf"), 
                width = 14, height = 7)
            print(f)
            dev.off()
        })
        
        try({
            plot <- VlnPlot(macrophages.integrated, features = gene, 
                group.by = "renamed.macrophages.subclusters", 
                split.by = "sample.id", pt.size = 0, combine = FALSE, 
                cols = colors.samples)
            pdf(paste0(outdir, "/27_", gene, "_VlnPlot_group_by_cluster_split_by_sample.pdf"), 
                width = 11, height = 8.5)
            print(plot)
            dev.off()
        })
    })
}
## [[1]]
## [[1]]
## [[1]]
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
try({
    d <- DotPlot(macrophages.integrated, features = goi, cols = c("#03539C", 
        "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.subclusters.sample.id") + 
        RotatedAxis()
    pdf(paste0(outdir, "/28_", paste(goi, sep = "_"), "_DotPlot_groupby_annotation.1.sample.id.pdf"), 
        width = (4 + length(goi) * 15/50), height = 17)
    print(d)
    dev.off()
})
## png 
##   2

Plots for more genes of interest: F480 (antibody), CD11b, CD11c, CD64, MHCII

In mouse: F4/80 (antibody) > Emr1, Adgre1 CD11b -> Itgam CD11c -> Itgax CD64 -> Fcgr1 MHCII -> H2-Ab

goi <- c("Adgre1", "Itgam", "Itgax", "Fcgr1", "H2-Ab1")
DefaultAssay(macrophages.integrated) <- "RNA"
for (gene in goi) {
    c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.", 
        node_colour = gene, node_colour_aggr = "mean")
    pdf(paste0(outdir, "/29_clustree_", gene, "_macrophages.pdf"))
    print(c)
    dev.off()
}
for (gene in goi) {
    try({
        Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
        try({
            f <- FeaturePlot(macrophages.integrated, features = gene, 
                split.by = "sample.id", cols = c("grey", "#E82564"), 
                order = T, min.cutoff = "q10", max.cutoff = "q90", 
                label = T)
            pdf(paste0(outdir, "/30_", gene, "_FeaturePlot_split_by_sample.pdf"), 
                width = 14, height = 7)
            print(f)
            dev.off()
        })
        
        try({
            plot <- VlnPlot(macrophages.integrated, features = gene, 
                group.by = "renamed.macrophages.subclusters", 
                split.by = "sample.id", pt.size = 0, combine = FALSE, 
                cols = colors.samples)
            pdf(paste0(outdir, "/31_", gene, "_VlnPlot_group_by_cluster_split_by_sample.pdf"), 
                width = 11, height = 8.5)
            print(plot)
            dev.off()
        })
    })
}
## [[1]]
## [[1]]
## [[1]]
## [[1]]
## [[1]]
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
try({
    d <- DotPlot(macrophages.integrated, features = goi, cols = c("#03539C", 
        "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.subclusters.sample.id") + 
        RotatedAxis()
    pdf(paste0(outdir, "/32_", paste(goi, sep = "_"), "_DotPlot_groupby_macrophages_subclusters.sample.id.pdf"), 
        width = (4 + length(goi) * 15/50), height = 17)
    print(d)
    dev.off()
})
## png 
##   2

Expression of marker Tim3/Havcr2 from reference Hasegawa 2021

https://jasn.asnjournals.org/content/early/2021/04/17/ASN.2020121723

DefaultAssay(macrophages.integrated) <- "RNA"
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
f <- FeaturePlot(macrophages.integrated, features = "Havcr2", 
    min.cutoff = "q9", order = T, label = T, label.size = 8, 
    repel = T)
f

pdf(paste0(outdir, "/33_FeaturePlot_Havcr2_macrophages.pdf"), 
    width = 7, height = 7)
f
dev.off()
## png 
##   2

QC metrics to possibly remove high mitochondrial expression cells

v <- VlnPlot(macrophages.integrated, features = c("nFeature_RNA", 
    "nCount_RNA", "percent.mito"), ncol = 3, cols = colors.samples, 
    group.by = "sample.id")
v

pdf(paste0(outdir, "/34_Seurat_QC_VlnPlot_samples.pdf"), width = 21, 
    height = 7)
v
dev.off()
## png 
##   2
summary(macrophages.integrated$percent.mito)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.7848  4.7934  7.9638 12.3109 16.1364 54.6878

20% seems like a reasonable threshold Jamie: Regarding your other question, I do think having a lower threshold for mitochondrial genes would be good in the macrophage subpopulations. I think the high mitochondrial content is mainly applicable to the tubular epithelial clusters.

First clean up of the subset by removing:

- genes that are not expressed in any cell

dim(macrophages.integrated)
## [1] 22399   588
keep_feature <- rowSums(as.matrix(GetAssayData(macrophages.integrated, 
    slot = "counts"))) > 0
summary(keep_feature)
##    Mode   FALSE    TRUE 
## logical    8106   14293

- cells in which no gene expression is detected

keep_cell <- colSums(as.matrix(GetAssayData(macrophages.integrated, 
    slot = "counts"))) > 0
summary(keep_cell)
##    Mode    TRUE 
## logical     588
macrophages.integrated.filtered <- macrophages.integrated[keep_feature, 
    keep_cell]
dim(macrophages.integrated.filtered)
## [1] 14293   588

Remove cells with mito expn > 20%

dim(macrophages.integrated.filtered)
## [1] 14293   588
keep_cell <- macrophages.integrated.filtered$percent.mito <= 
    20
summary(keep_cell)
##    Mode   FALSE    TRUE 
## logical     108     480
macrophages.integrated.filtered <- macrophages.integrated.filtered[, 
    keep_cell]
dim(macrophages.integrated.filtered)
## [1] 14293   480

New subclsutering after filtering

DefaultAssay(macrophages.integrated.filtered) <- "RNA"
macrophages.2.list <- SplitObject(macrophages.integrated.filtered, 
    split.by = "sample.id")
macrophages.2.list
## $C0
## An object of class Seurat 
## 16293 features across 212 samples within 2 assays 
## Active assay: RNA (14293 features, 0 variable features)
##  1 other assay present: integrated
##  2 dimensional reductions calculated: pca, umap
## 
## $C24
## An object of class Seurat 
## 16293 features across 268 samples within 2 assays 
## Active assay: RNA (14293 features, 0 variable features)
##  1 other assay present: integrated
##  2 dimensional reductions calculated: pca, umap
for (i in 1:length(macrophages.2.list)) {
    macrophages.2.list[[i]] <- FindVariableFeatures(macrophages.2.list[[i]], 
        selection.method = "vst", nfeatures = 2000, verbose = FALSE)
}
k.filter <- min(200, min(sapply(macrophages.2.list, ncol)))
k.filter
## [1] 200
macrophages.2.anchors <- FindIntegrationAnchors(object.list = macrophages.2.list, 
    dims = 1:30, k.filter = k.filter)
macrophages.2.integrated <- IntegrateData(anchorset = macrophages.2.anchors, 
    dims = 1:30)
DefaultAssay(macrophages.2.integrated) <- "integrated"
macrophages.2.integrated <- ScaleData(macrophages.2.integrated, 
    verbose = T)
dim(macrophages.2.integrated)
## [1] 2000  480
macrophages.2.integrated <- RunPCA(macrophages.2.integrated, 
    npcs = 50, verbose = T)
system.time(macrophages.2.integrated <- JackStraw(macrophages.2.integrated, 
    num.replicate = 100, dims = 50))
##    user  system elapsed 
##  63.442   0.093  63.595
macrophages.2.integrated <- ScoreJackStraw(macrophages.2.integrated, 
    dims = 1:50)
j <- JackStrawPlot(macrophages.2.integrated, dims = 1:50)
pdf(paste0(outdir, "/35_JackStrawPlot_macrophages_2.pdf"))
j
dev.off()
## png 
##   2
macrophages.2.integrated <- FindNeighbors(macrophages.2.integrated, 
    dims = 1:48)
macrophages.2.integrated <- RunUMAP(macrophages.2.integrated, 
    dims = 1:48)
for (i in seq(0, 2, 0.1)) {
    macrophages.2.integrated <- FindClusters(macrophages.2.integrated, 
        resolution = i, verbose = F)
    macrophages.2.integrated[[paste0("macrophages_2_subclusters_res.", 
        i)]] <- macrophages.2.integrated$seurat_clusters
}
# install.packages('clustree')
library(clustree)
c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.")
c

pdf(paste0(outdir, "/36_clustree_macrophages_2.pdf"))
c
dev.off()
## png 
##   2

Using the stability index to assess clusters The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions and is automatically calculated when a clustering tree is built.

c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.", 
    node_colour = "sc3_stability")
pdf(paste0(outdir, "/37_clustree_macrophages_2_sc3_stability.pdf"))
c
dev.off()
## png 
##   2

SOCS3, IL-1R2, IL1rn

Socs3, Il1r2, Il1rn

c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.", 
    node_colour = "Socs3", node_colour_aggr = "mean")
c

pdf(paste0(outdir, "/38_clustree_macrophages_2_Socs3.pdf"))
c
dev.off()
## png 
##   2
c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.", 
    node_colour = "Il1r2", node_colour_aggr = "mean")
c

pdf(paste0(outdir, "/39_clustree_macrophages_2_Il1r2.pdf"))
c
dev.off()
## png 
##   2
c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.", 
    node_colour = "Il1rn", node_colour_aggr = "mean")
c

pdf(paste0(outdir, "/40_clustree_macrophages_2_Il1rn.pdf"))
c
dev.off()
## png 
##   2

RES 0.2

# install.packages('Polychrome')
library(Polychrome)
set.seed(5658807)  # for reproducibility
levels <- levels(as.factor(macrophages.2.integrated$macrophages_2_subclusters_res.0.2))
rainbow2.6c <- c("#03539C", "#12999E", "#B7CE05", "#FAEB09", 
    "#FDA908", "#E82564")
colors.macrophages.2.subclusters <- createPalette(length(levels), 
    rainbow2.6c, M = 1000)
names(colors.macrophages.2.subclusters) <- levels
colors.macrophages.2.subclusters
##         0         1         2      <NA>      <NA>      <NA>      <NA>      <NA> 
## "#1662B5" "#0D9A9D" "#C2DB1C" "#FBE22E" "#FCA700" "#F33877" "#FA35F5" "#FFB1DB" 
##      <NA>      <NA>      <NA> 
## "#40F897" "#63691C" "#4F00FF"
slices <- rep(1, length(colors.macrophages.2.subclusters))
pie(slices, col = colors.macrophages.2.subclusters, labels = names(colors.macrophages.2.subclusters))

Idents(macrophages.2.integrated) <- "macrophages_2_subclusters_res.0.2"
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 8, cols = colors.macrophages.2.subclusters)

pdf(paste0(outdir, "/41_DimPlot_umap_macrophages_2_subclusters_res.0.2.pdf"))
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 6, cols = colors.macrophages.2.subclusters)
dev.off()
## png 
##   2
colors.samples <- c("#12999E", "#FDA908")
names(colors.samples) <- c("C0", "C24")
macrophages.2.integrated$sample.id <- factor(macrophages.2.integrated$sample.id, 
    levels = names(colors.samples))
p <- DimPlot(macrophages.2.integrated, reduction = "umap", group.by = "sample.id", 
    cols = colors.samples, label = F)
pdf(paste0(outdir, "/42_UMAP_macrophages_2_samples.pdf"), width = 10, 
    height = 9)
p
dev.off()
## png 
##   2
pdf(paste0(outdir, "/43_UMAP_macrophages_2_samples_noLegend.pdf"), 
    width = 7, height = 7)
p + NoLegend()
dev.off()
## png 
##   2

Rename subclusters

new.cluster.ids <- paste0("Macro", 1:length(levels(as.factor(macrophages.2.integrated$macrophages_2_subclusters_res.0.2))))
names(new.cluster.ids) <- levels(macrophages.2.integrated)
macrophages.2.integrated <- RenameIdents(macrophages.2.integrated, 
    new.cluster.ids)
macrophages.2.integrated[["renamed.macrophages.2.subclusters"]] <- macrophages.2.integrated@active.ident
Idents(macrophages.2.integrated) <- "renamed.macrophages.2.subclusters"
names(colors.macrophages.2.subclusters) <- levels(macrophages.2.integrated$renamed.macrophages.2.subclusters)
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 8, cols = colors.macrophages.2.subclusters)

pdf(paste0(outdir, "/44_DimPlot_umap_renamed_macrophages_2_subclusters.pdf"))
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 6, cols = colors.macrophages.2.subclusters)
dev.off()
## png 
##   2
Idents(macrophages.2.integrated) <- "renamed.macrophages.2.subclusters"
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1, 
    label = F, label.size = 8, cols = colors.macrophages.2.subclusters) + 
    NoLegend()

pdf(paste0(outdir, "/45_DimPlot_umap_renamed_macrophages_2_subclusters_noLabels_noLegend.pdf"))
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1, 
    label = F, label.size = 6, cols = colors.macrophages.2.subclusters) + 
    NoLegend()
dev.off()
## png 
##   2
saveRDS(macrophages.2.integrated, paste0(outdir, "/46.macrophages.2.integrated.rds"))

Conserved markers

t <- table(macrophages.2.integrated$renamed.macrophages.2.subclusters, 
    macrophages.2.integrated$sample.id)
t
##         
##           C0 C24
##   Macro1 139 151
##   Macro2  38  75
##   Macro3  35  42
openxlsx::write.xlsx(as.data.frame.matrix(t), file = paste0(outdir, 
    "/47_Sample_macrophages_2_subcluster_composition.xlsx"), 
    rowNames = T, colNames = T)
DefaultAssay(macrophages.2.integrated) <- "RNA"
Idents(macrophages.2.integrated) <- "renamed.macrophages.2.subclusters"
subclusters <- levels(as.factor(macrophages.2.integrated$renamed.macrophages.2.subclusters))
conserved.markers <- data.frame(matrix(ncol = 14))
for (c in subclusters) {
    print(c)
    markers.c <- FindConservedMarkers(macrophages.2.integrated, 
        ident.1 = c, grouping.var = "sample.id", verbose = T, 
        logfc.threshold = -Inf, min.pct = -Inf, min.cells.feature = -Inf, 
        min.cells.group = -Inf)
    markers.c <- cbind(data.frame(cluster = rep(c, dim(markers.c)[1]), 
        gene = rownames(markers.c)), markers.c)
    write.table(markers.c, file = paste0(outdir, "/48_markers_", 
        c, "_macrophages_2_subclusters.txt"))
    colnames(conserved.markers) <- colnames(markers.c)
    conserved.markers <- rbind(conserved.markers, markers.c)
}
## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"
# for (c in 1:4) { subcluster <- paste0('AM', c) markers.c <-
# read.table(file = paste0(outdir, '/31_markers_',
# subcluster, '_AM_subclusters.txt'))
# colnames(conserved.markers) <- colnames(markers.c)
# conserved.markers <- rbind(conserved.markers, markers.c) }
# head(conserved.markers)
# markers.AM5 <- read.table(file = paste0(outdir,
# '/31_markers_', 'AM5', '_AM_subclusters.txt'))
# colnames(markers.AM5) colnames(conserved.markers)
# for ( i in 1:10 ) { markers.AM5[, 34+i] <- rep(NA,
# dim(markers.AM5)[1]) } markers.AM5 <- markers.AM5[, c(1:12,
# 35:44, 13:34)] colnames(markers.AM5) <-
# colnames(conserved.markers) head(markers.AM5)
# conserved.markers <- rbind(conserved.markers, markers.AM5)
conserved.markers <- conserved.markers[-1, ]
conserved.markers <- conserved.markers[, c(1:2, 13:14, 3:12)]
head(conserved.markers)
##        cluster   gene     max_pval minimump_p_val     C0_p_val C0_avg_log2FC
## Sdc4    Macro1   Sdc4 2.202498e-01   9.940550e-25 2.202498e-01   -0.07585808
## Cd14    Macro1   Cd14 3.933950e-13   2.413485e-23 3.933950e-13    1.97725880
## Cebpb   Macro1  Cebpb 1.065972e-02   1.624807e-20 1.065972e-02    0.66786332
## Pla2g7  Macro1 Pla2g7 1.069071e-05   1.458113e-18 1.069071e-05    1.31291835
## Pid1    Macro1   Pid1 1.655140e-01   4.899195e-18 1.655140e-01    0.04600540
## Ctsb    Macro1   Ctsb 1.454003e-06   1.396807e-17 1.454003e-06    0.84584724
##        C0_pct.1 C0_pct.2 C0_p_val_adj    C24_p_val C24_avg_log2FC C24_pct.1
## Sdc4      0.309    0.219 1.000000e+00 4.970275e-25       1.494107     1.000
## Cd14      0.755    0.233 5.622794e-09 1.206743e-23       1.878361     0.980
## Cebpb     0.194    0.068 1.000000e+00 8.124037e-21       1.492180     0.987
## Pla2g7    0.324    0.055 1.528024e-01 7.290567e-19       1.376837     0.709
## Pid1      0.460    0.384 1.000000e+00 2.449597e-18       1.276004     0.788
## Ctsb      0.806    0.534 2.078207e-02 6.984035e-18       1.269536     0.967
##        C24_pct.2 C24_p_val_adj
## Sdc4       0.667  7.104014e-21
## Cd14       0.496  1.724797e-19
## Cebpb      0.556  1.161169e-16
## Pla2g7     0.103  1.042041e-14
## Pid1       0.231  3.501209e-14
## Ctsb       0.556  9.982281e-14

Only top markers that are positive markers

conserved.markers$marker.type <- apply(conserved.markers, 1, function(x) {
  y <- as.numeric(x)
  if ( (if (is.na(y[6])) {TRUE} else {y[6]>0})
       & (if (is.na(y[11])) {TRUE} else {y[11]>0})
       # & (if (is.na(y[16])) {TRUE} else {y[16]>0})
       # & (if (is.na(y[21])) {TRUE} else {y[21]>0})
       # & (if (is.na(y[26])) {TRUE} else {y[26]>0})
       # & (if (is.na(y[31])) {TRUE} else {y[31]>0})
       # & (if (is.na(y[36])) {TRUE} else {y[36]>0})
       # & (if (is.na(y[41])) {TRUE} else {y[41]>0})
       )
    {"POS"}
  else if ( ( if (is.na(y[6])) {TRUE} else {y[6]<0})
       & (if (is.na(y[11])) {TRUE} else {y[11]<0})
       # & (if (is.na(y[16])) {TRUE} else {y[16]<0})
       # & (if (is.na(y[21])) {TRUE} else {y[21]<0})
       # & (if (is.na(y[26])) {TRUE} else {y[26]<0})
       # & (if (is.na(y[31])) {TRUE} else {y[31]<0})
       # & (if (is.na(y[36])) {TRUE} else {y[36]<0})
       # & (if (is.na(y[41])) {TRUE} else {y[41]<0})
       ) 
    {"NEG"}
  else {"UND"}
  })
conserved.markers <- conserved.markers[, c(1:4, 15, 5:14)]
openxlsx::write.xlsx(conserved.markers, paste0(outdir, "/49_conserved_markers_macrophages_2.xlsx"), 
    colNames = T)
head(conserved.markers)
##        cluster   gene     max_pval minimump_p_val marker.type     C0_p_val
## Sdc4    Macro1   Sdc4 2.202498e-01   9.940550e-25         UND 2.202498e-01
## Cd14    Macro1   Cd14 3.933950e-13   2.413485e-23         POS 3.933950e-13
## Cebpb   Macro1  Cebpb 1.065972e-02   1.624807e-20         POS 1.065972e-02
## Pla2g7  Macro1 Pla2g7 1.069071e-05   1.458113e-18         POS 1.069071e-05
## Pid1    Macro1   Pid1 1.655140e-01   4.899195e-18         POS 1.655140e-01
## Ctsb    Macro1   Ctsb 1.454003e-06   1.396807e-17         POS 1.454003e-06
##        C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj    C24_p_val C24_avg_log2FC
## Sdc4     -0.07585808    0.309    0.219 1.000000e+00 4.970275e-25       1.494107
## Cd14      1.97725880    0.755    0.233 5.622794e-09 1.206743e-23       1.878361
## Cebpb     0.66786332    0.194    0.068 1.000000e+00 8.124037e-21       1.492180
## Pla2g7    1.31291835    0.324    0.055 1.528024e-01 7.290567e-19       1.376837
## Pid1      0.04600540    0.460    0.384 1.000000e+00 2.449597e-18       1.276004
## Ctsb      0.84584724    0.806    0.534 2.078207e-02 6.984035e-18       1.269536
##        C24_pct.1 C24_pct.2 C24_p_val_adj
## Sdc4       1.000     0.667  7.104014e-21
## Cd14       0.980     0.496  1.724797e-19
## Cebpb      0.987     0.556  1.161169e-16
## Pla2g7     0.709     0.103  1.042041e-14
## Pid1       0.788     0.231  3.501209e-14
## Ctsb       0.967     0.556  9.982281e-14
saveRDS(macrophages.2.integrated, paste0(outdir, "/50.macrophages.2.integrated.rds"))
conserved.markers.pos <- conserved.markers[conserved.markers$marker.type == 
    "POS", ]
head(conserved.markers.pos)
##        cluster   gene     max_pval minimump_p_val marker.type     C0_p_val
## Cd14    Macro1   Cd14 3.933950e-13   2.413485e-23         POS 3.933950e-13
## Cebpb   Macro1  Cebpb 1.065972e-02   1.624807e-20         POS 1.065972e-02
## Pla2g7  Macro1 Pla2g7 1.069071e-05   1.458113e-18         POS 1.069071e-05
## Pid1    Macro1   Pid1 1.655140e-01   4.899195e-18         POS 1.655140e-01
## Ctsb    Macro1   Ctsb 1.454003e-06   1.396807e-17         POS 1.454003e-06
## C1qa    Macro1   C1qa 5.312464e-13   1.480332e-16         POS 7.401660e-17
##        C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj    C24_p_val C24_avg_log2FC
## Cd14       1.9772588    0.755    0.233 5.622794e-09 1.206743e-23       1.878361
## Cebpb      0.6678633    0.194    0.068 1.000000e+00 8.124037e-21       1.492180
## Pla2g7     1.3129183    0.324    0.055 1.528024e-01 7.290567e-19       1.376837
## Pid1       0.0460054    0.460    0.384 1.000000e+00 2.449597e-18       1.276004
## Ctsb       0.8458472    0.806    0.534 2.078207e-02 6.984035e-18       1.269536
## C1qa       1.6955557    0.957    0.452 1.057919e-12 5.312464e-13       1.032159
##        C24_pct.1 C24_pct.2 C24_p_val_adj
## Cd14       0.980     0.496  1.724797e-19
## Cebpb      0.987     0.556  1.161169e-16
## Pla2g7     0.709     0.103  1.042041e-14
## Pid1       0.788     0.231  3.501209e-14
## Ctsb       0.967     0.556  9.982281e-14
## C1qa       0.967     0.538  7.593105e-09
table(macrophages.2.integrated$sample.id)
## 
##  C0 C24 
## 212 268
colnames(conserved.markers.pos)
##  [1] "cluster"        "gene"           "max_pval"       "minimump_p_val"
##  [5] "marker.type"    "C0_p_val"       "C0_avg_log2FC"  "C0_pct.1"      
##  [9] "C0_pct.2"       "C0_p_val_adj"   "C24_p_val"      "C24_avg_log2FC"
## [13] "C24_pct.1"      "C24_pct.2"      "C24_p_val_adj"
library(dplyr)
top1.pos <- conserved.markers.pos %>% group_by(cluster) %>% top_n(n = 1, 
    wt = -minimump_p_val)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = -max_pval)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = C24_avg_log2FC)
head(top1.pos)
## # A tibble: 3 x 15
## # Groups:   cluster [3]
##   cluster gene  max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
##   <chr>   <chr>    <dbl>          <dbl> <chr>          <dbl>         <dbl>
## 1 Macro1  Cd14  3.93e-13       2.41e-23 POS         3.93e-13         1.98 
## 2 Macro2  Gpx3  7.94e- 8       3.02e-21 POS         7.94e- 8         1.24 
## 3 Macro3  Pfkp  3.18e- 2       1.29e-31 POS         3.18e- 2         0.168
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## #   C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## #   C24_p_val_adj <dbl>
DefaultAssay(macrophages.2.integrated) <- "RNA"
f <- FeaturePlot(macrophages.2.integrated, features = top1.pos$gene, 
    min.cutoff = "q9", order = T, label = T, label.size = 8, 
    repel = T)
f

pdf(paste0(outdir, "/51_FeaturePlot_top_pos_markers_macrophages_2.pdf"), 
    width = 14, height = 14)
f
dev.off()
## png 
##   2

DotPlot of top markers

DefaultAssay(macrophages.2.integrated) <- "RNA"
d <- DotPlot(object = macrophages.2.integrated, features = top1.pos$gene, 
    cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.2.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/52_DotPlot_top_cell_types_pos_markers_clusters.macrophages_2.pdf"), 
    width = 7, height = 7)
d
dev.off()
## png 
##   2

DotPlot of top 10 markers

library(dplyr)
top10.pos <- conserved.markers.pos %>% group_by(cluster) %>% 
    top_n(n = 10, wt = -minimump_p_val)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10, 
    wt = -max_pval)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10, 
    wt = C24_avg_log2FC)
head(top10.pos)
## # A tibble: 6 x 15
## # Groups:   cluster [1]
##   cluster gene  max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
##   <chr>   <chr>    <dbl>          <dbl> <chr>          <dbl>         <dbl>
## 1 Macro1  Cd14  3.93e-13       2.41e-23 POS         3.93e-13        1.98  
## 2 Macro1  Cebpb 1.07e- 2       1.62e-20 POS         1.07e- 2        0.668 
## 3 Macro1  Pla2… 1.07e- 5       1.46e-18 POS         1.07e- 5        1.31  
## 4 Macro1  Pid1  1.66e- 1       4.90e-18 POS         1.66e- 1        0.0460
## 5 Macro1  Ctsb  1.45e- 6       1.40e-17 POS         1.45e- 6        0.846 
## 6 Macro1  C1qa  5.31e-13       1.48e-16 POS         7.40e-17        1.70  
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## #   C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## #   C24_p_val_adj <dbl>
DefaultAssay(macrophages.2.integrated) <- "RNA"
d <- DotPlot(object = macrophages.2.integrated, features = top10.pos$gene, 
    cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.2.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/53_DotPlot_top10_cell_types_pos_markers_clusters.macrophages.2.pdf"), 
    width = 14, height = 7)
d
dev.off()
## png 
##   2

Segregation of clusters by various sources of uninteresting variation

Explore additional metrics, such as the number of UMIs and genes per cell, and mitochondrial gene expression by UMAP.

# Determine metrics to plot present in
# tcells.integrated@meta.data
metrics <- c("nCount_RNA", "nFeature_RNA", "percent.mito")

f <- FeaturePlot(macrophages.2.integrated, reduction = "umap", 
    features = metrics, pt.size = 0.4, order = TRUE, min.cutoff = "q10", 
    label = TRUE)
f

pdf(paste0(outdir, "/54_FeaturePlot_macrophages.2.integrated_umap_metrics.pdf"), 
    width = 14, height = 14)
f
dev.off()
## png 
##   2

Remove cells with mito expn > 10%

dim(macrophages.integrated.filtered)
## [1] 14293   480
keep_cell <- macrophages.integrated.filtered$percent.mito <= 
    10
summary(keep_cell)
##    Mode   FALSE    TRUE 
## logical     132     348
macrophages.integrated.filtered <- macrophages.integrated.filtered[, 
    keep_cell]
dim(macrophages.integrated.filtered)
## [1] 14293   348

New subclustering after filtering

DefaultAssay(macrophages.integrated.filtered) <- "RNA"
macrophages.3.list <- SplitObject(macrophages.integrated.filtered, 
    split.by = "sample.id")
macrophages.3.list
## $C0
## An object of class Seurat 
## 16293 features across 156 samples within 2 assays 
## Active assay: RNA (14293 features, 0 variable features)
##  1 other assay present: integrated
##  2 dimensional reductions calculated: pca, umap
## 
## $C24
## An object of class Seurat 
## 16293 features across 192 samples within 2 assays 
## Active assay: RNA (14293 features, 0 variable features)
##  1 other assay present: integrated
##  2 dimensional reductions calculated: pca, umap
for (i in 1:length(macrophages.3.list)) {
    macrophages.3.list[[i]] <- FindVariableFeatures(macrophages.3.list[[i]], 
        selection.method = "vst", nfeatures = 2000, verbose = FALSE)
}
k.filter <- min(200, min(sapply(macrophages.3.list, ncol)))
k.filter
## [1] 156
macrophages.3.anchors <- FindIntegrationAnchors(object.list = macrophages.3.list, 
    dims = 1:30, k.filter = k.filter)
macrophages.3.integrated <- IntegrateData(anchorset = macrophages.3.anchors, 
    dims = 1:30)
DefaultAssay(macrophages.3.integrated) <- "integrated"
macrophages.3.integrated <- ScaleData(macrophages.3.integrated, 
    verbose = T)
dim(macrophages.3.integrated)
## [1] 2000  348
macrophages.3.integrated <- RunPCA(macrophages.3.integrated, 
    npcs = 50, verbose = T)
system.time(macrophages.3.integrated <- JackStraw(macrophages.3.integrated, 
    num.replicate = 100, dims = 50))
##    user  system elapsed 
##  57.091   0.007  57.147
macrophages.3.integrated <- ScoreJackStraw(macrophages.3.integrated, 
    dims = 1:50)
j <- JackStrawPlot(macrophages.3.integrated, dims = 1:50)
pdf(paste0(outdir, "/55_JackStrawPlot_macrophages_3.pdf"))
j
dev.off()
## png 
##   2
macrophages.3.integrated <- FindNeighbors(macrophages.3.integrated, 
    dims = 1:42)
macrophages.3.integrated <- RunUMAP(macrophages.3.integrated, 
    dims = 1:42)
for (i in seq(0, 2, 0.1)) {
    macrophages.3.integrated <- FindClusters(macrophages.3.integrated, 
        resolution = i, verbose = F)
    macrophages.3.integrated[[paste0("macrophages_3_subclusters_res.", 
        i)]] <- macrophages.3.integrated$seurat_clusters
}
# install.packages('clustree')
library(clustree)
c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.")
c

pdf(paste0(outdir, "/56_clustree_macrophages_3.pdf"))
c
dev.off()
## png 
##   2

Using the stability index to assess clusters The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions and is automatically calculated when a clustering tree is built.

c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.", 
    node_colour = "sc3_stability")
pdf(paste0(outdir, "/57_clustree_macrophages_3_sc3_stability.pdf"))
c
dev.off()
## png 
##   2

SOCS3, IL-1R2, IL1rn

Socs3, Il1r2, Il1rn

c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.", 
    node_colour = "Socs3", node_colour_aggr = "mean")
c

pdf(paste0(outdir, "/58_clustree_macrophages_3_Socs3.pdf"))
c
dev.off()
## png 
##   2
c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.", 
    node_colour = "Il1r2", node_colour_aggr = "mean")
c

pdf(paste0(outdir, "/59_clustree_macrophages_3_Il1r2.pdf"))
c
dev.off()
## png 
##   2
c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.", 
    node_colour = "Il1rn", node_colour_aggr = "mean")
c

pdf(paste0(outdir, "/60_clustree_macrophages_3_Il1rn.pdf"))
c
dev.off()
## png 
##   2

RES 0.3

# install.packages('Polychrome')
library(Polychrome)
set.seed(5658807)  # for reproducibility
levels <- levels(as.factor(macrophages.3.integrated$macrophages_3_subclusters_res.0.3))
rainbow2.6c <- c("#03539C", "#12999E", "#B7CE05", "#FAEB09", 
    "#FDA908", "#E82564")
colors.macrophages.3.subclusters <- createPalette(length(levels), 
    rainbow2.6c, M = 1000)
names(colors.macrophages.3.subclusters) <- levels
colors.macrophages.3.subclusters
##         0         1         2      <NA>      <NA>      <NA>      <NA>      <NA> 
## "#1662B5" "#0D9A9D" "#C2DB1C" "#FBE22E" "#FCA700" "#F33877" "#FA35F5" "#FFB1DB" 
##      <NA>      <NA>      <NA> 
## "#40F897" "#63691C" "#4F00FF"
slices <- rep(1, length(colors.macrophages.3.subclusters))
pie(slices, col = colors.macrophages.3.subclusters, labels = names(colors.macrophages.3.subclusters))

Idents(macrophages.3.integrated) <- "macrophages_3_subclusters_res.0.3"
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 8, cols = colors.macrophages.3.subclusters)

pdf(paste0(outdir, "/61_DimPlot_umap_macrophages_3_subclusters_res.0.3.pdf"))
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 6, cols = colors.macrophages.3.subclusters)
dev.off()
## png 
##   2
colors.samples <- c("#12999E", "#FDA908")
names(colors.samples) <- c("C0", "C24")
macrophages.3.integrated$sample.id <- factor(macrophages.3.integrated$sample.id, 
    levels = names(colors.samples))
p <- DimPlot(macrophages.3.integrated, reduction = "umap", group.by = "sample.id", 
    cols = colors.samples, label = F)
pdf(paste0(outdir, "/62_UMAP_macrophages_3_samples.pdf"), width = 10, 
    height = 9)
p
dev.off()
## png 
##   2
pdf(paste0(outdir, "/63_UMAP_macrophages_3_samples_noLegend.pdf"), 
    width = 7, height = 7)
p + NoLegend()
dev.off()
## png 
##   2

Rename subclusters

new.cluster.ids <- paste0("Macro", 1:length(levels(as.factor(macrophages.3.integrated$macrophages_3_subclusters_res.0.3))))
names(new.cluster.ids) <- levels(macrophages.3.integrated)
macrophages.3.integrated <- RenameIdents(macrophages.3.integrated, 
    new.cluster.ids)
macrophages.3.integrated[["renamed.macrophages.3.subclusters"]] <- macrophages.3.integrated@active.ident
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
names(colors.macrophages.3.subclusters) <- levels(macrophages.3.integrated$renamed.macrophages.3.subclusters)
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 8, cols = colors.macrophages.3.subclusters)

pdf(paste0(outdir, "/64_DimPlot_umap_renamed_macrophages_3_subclusters.pdf"))
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1, 
    label = T, label.size = 6, cols = colors.macrophages.3.subclusters)
dev.off()
## png 
##   2
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1, 
    label = F, label.size = 8, cols = colors.macrophages.3.subclusters) + 
    NoLegend()

pdf(paste0(outdir, "/65_DimPlot_umap_renamed_macrophages_3_subclusters_noLabels_noLegend.pdf"))
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1, 
    label = F, label.size = 6, cols = colors.macrophages.3.subclusters) + 
    NoLegend()
dev.off()
## png 
##   2
saveRDS(macrophages.3.integrated, paste0(outdir, "/66.macrophages.3.integrated.rds"))

Conserved markers

t <- table(macrophages.3.integrated$renamed.macrophages.3.subclusters, 
    macrophages.3.integrated$sample.id)
t
##         
##           C0 C24
##   Macro1 108 118
##   Macro2  34  42
##   Macro3  14  32
openxlsx::write.xlsx(as.data.frame.matrix(t), file = paste0(outdir, 
    "/67_Sample_macrophages_3_subcluster_composition.xlsx"), 
    rowNames = T, colNames = T)
DefaultAssay(macrophages.3.integrated) <- "RNA"
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
subclusters <- levels(as.factor(macrophages.3.integrated$renamed.macrophages.3.subclusters))
conserved.markers <- data.frame(matrix(ncol = 14))
for (c in subclusters) {
    print(c)
    markers.c <- FindConservedMarkers(macrophages.3.integrated, 
        ident.1 = c, grouping.var = "sample.id", verbose = T, 
        logfc.threshold = -Inf, min.pct = -Inf, min.cells.feature = -Inf, 
        min.cells.group = -Inf)
    markers.c <- cbind(data.frame(cluster = rep(c, dim(markers.c)[1]), 
        gene = rownames(markers.c)), markers.c)
    write.table(markers.c, file = paste0(outdir, "/68_markers_", 
        c, "_macrophages_3_subclusters.txt"))
    colnames(conserved.markers) <- colnames(markers.c)
    conserved.markers <- rbind(conserved.markers, markers.c)
}
## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"
# for (c in 1:4) { subcluster <- paste0('AM', c) markers.c <-
# read.table(file = paste0(outdir, '/31_markers_',
# subcluster, '_AM_subclusters.txt'))
# colnames(conserved.markers) <- colnames(markers.c)
# conserved.markers <- rbind(conserved.markers, markers.c) }
# head(conserved.markers)
# markers.AM5 <- read.table(file = paste0(outdir,
# '/31_markers_', 'AM5', '_AM_subclusters.txt'))
# colnames(markers.AM5) colnames(conserved.markers)
# for ( i in 1:10 ) { markers.AM5[, 34+i] <- rep(NA,
# dim(markers.AM5)[1]) } markers.AM5 <- markers.AM5[, c(1:12,
# 35:44, 13:34)] colnames(markers.AM5) <-
# colnames(conserved.markers) head(markers.AM5)
# conserved.markers <- rbind(conserved.markers, markers.AM5)
conserved.markers <- conserved.markers[-1, ]
conserved.markers <- conserved.markers[, c(1:2, 13:14, 3:12)]
head(conserved.markers)
##       cluster  gene     max_pval minimump_p_val     C0_p_val C0_avg_log2FC
## Cd14   Macro1  Cd14 3.903045e-10   2.835224e-21 3.903045e-10     1.9896992
## C1qa   Macro1  C1qa 7.937402e-15   9.276577e-17 4.638289e-17     2.2670625
## C1qc   Macro1  C1qc 7.233602e-13   3.884857e-15 7.233602e-13     1.6955969
## Cebpb  Macro1 Cebpb 3.638433e-01   7.452346e-15 3.638433e-01     0.1487639
## Sdc4   Macro1  Sdc4 4.119571e-02   1.039183e-14 4.119571e-02     0.3247144
## C1qb   Macro1  C1qb 8.047736e-12   2.132823e-14 1.066411e-14     1.9546166
##       C0_pct.1 C0_pct.2 C0_p_val_adj    C24_p_val C24_avg_log2FC C24_pct.1
## Cd14     0.796    0.271 5.578622e-06 1.417612e-21       2.201101     1.000
## C1qa     0.972    0.396 6.629506e-13 7.937402e-15       1.509209     0.992
## C1qc     0.981    0.375 1.033899e-08 1.942429e-15       1.617483     0.975
## Cebpb    0.185    0.125 1.000000e+00 3.726173e-15       1.545982     0.983
## Sdc4     0.306    0.146 1.000000e+00 5.195913e-15       1.213547     1.000
## C1qb     0.972    0.417 1.524222e-10 8.047736e-12       1.113736     1.000
##       C24_pct.2 C24_p_val_adj
## Cd14      0.541  2.026193e-17
## C1qa      0.405  1.134493e-10
## C1qc      0.365  2.776313e-11
## Cebpb     0.581  5.325819e-11
## Sdc4      0.784  7.426519e-11
## C1qb      0.486  1.150263e-07

Only top markers that are positive markers

conserved.markers$marker.type <- apply(conserved.markers, 1, function(x) {
  y <- as.numeric(x)
  if ( (if (is.na(y[6])) {TRUE} else {y[6]>0})
       & (if (is.na(y[11])) {TRUE} else {y[11]>0})
       # & (if (is.na(y[16])) {TRUE} else {y[16]>0})
       # & (if (is.na(y[21])) {TRUE} else {y[21]>0})
       # & (if (is.na(y[26])) {TRUE} else {y[26]>0})
       # & (if (is.na(y[31])) {TRUE} else {y[31]>0})
       # & (if (is.na(y[36])) {TRUE} else {y[36]>0})
       # & (if (is.na(y[41])) {TRUE} else {y[41]>0})
       )
    {"POS"}
  else if ( ( if (is.na(y[6])) {TRUE} else {y[6]<0})
       & (if (is.na(y[11])) {TRUE} else {y[11]<0})
       # & (if (is.na(y[16])) {TRUE} else {y[16]<0})
       # & (if (is.na(y[21])) {TRUE} else {y[21]<0})
       # & (if (is.na(y[26])) {TRUE} else {y[26]<0})
       # & (if (is.na(y[31])) {TRUE} else {y[31]<0})
       # & (if (is.na(y[36])) {TRUE} else {y[36]<0})
       # & (if (is.na(y[41])) {TRUE} else {y[41]<0})
       ) 
    {"NEG"}
  else {"UND"}
  })
conserved.markers <- conserved.markers[, c(1:4, 15, 5:14)]
openxlsx::write.xlsx(conserved.markers, paste0(outdir, "/69_conserved_markers_macrophages_3.xlsx"), 
    colNames = T)
head(conserved.markers)
##       cluster  gene     max_pval minimump_p_val marker.type     C0_p_val
## Cd14   Macro1  Cd14 3.903045e-10   2.835224e-21         POS 3.903045e-10
## C1qa   Macro1  C1qa 7.937402e-15   9.276577e-17         POS 4.638289e-17
## C1qc   Macro1  C1qc 7.233602e-13   3.884857e-15         POS 7.233602e-13
## Cebpb  Macro1 Cebpb 3.638433e-01   7.452346e-15         POS 3.638433e-01
## Sdc4   Macro1  Sdc4 4.119571e-02   1.039183e-14         POS 4.119571e-02
## C1qb   Macro1  C1qb 8.047736e-12   2.132823e-14         POS 1.066411e-14
##       C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj    C24_p_val C24_avg_log2FC
## Cd14      1.9896992    0.796    0.271 5.578622e-06 1.417612e-21       2.201101
## C1qa      2.2670625    0.972    0.396 6.629506e-13 7.937402e-15       1.509209
## C1qc      1.6955969    0.981    0.375 1.033899e-08 1.942429e-15       1.617483
## Cebpb     0.1487639    0.185    0.125 1.000000e+00 3.726173e-15       1.545982
## Sdc4      0.3247144    0.306    0.146 1.000000e+00 5.195913e-15       1.213547
## C1qb      1.9546166    0.972    0.417 1.524222e-10 8.047736e-12       1.113736
##       C24_pct.1 C24_pct.2 C24_p_val_adj
## Cd14      1.000     0.541  2.026193e-17
## C1qa      0.992     0.405  1.134493e-10
## C1qc      0.975     0.365  2.776313e-11
## Cebpb     0.983     0.581  5.325819e-11
## Sdc4      1.000     0.784  7.426519e-11
## C1qb      1.000     0.486  1.150263e-07
saveRDS(macrophages.3.integrated, paste0(outdir, "/70.macrophages.3.integrated.rds"))
conserved.markers.pos <- conserved.markers[conserved.markers$marker.type == 
    "POS", ]
head(conserved.markers.pos)
##       cluster  gene     max_pval minimump_p_val marker.type     C0_p_val
## Cd14   Macro1  Cd14 3.903045e-10   2.835224e-21         POS 3.903045e-10
## C1qa   Macro1  C1qa 7.937402e-15   9.276577e-17         POS 4.638289e-17
## C1qc   Macro1  C1qc 7.233602e-13   3.884857e-15         POS 7.233602e-13
## Cebpb  Macro1 Cebpb 3.638433e-01   7.452346e-15         POS 3.638433e-01
## Sdc4   Macro1  Sdc4 4.119571e-02   1.039183e-14         POS 4.119571e-02
## C1qb   Macro1  C1qb 8.047736e-12   2.132823e-14         POS 1.066411e-14
##       C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj    C24_p_val C24_avg_log2FC
## Cd14      1.9896992    0.796    0.271 5.578622e-06 1.417612e-21       2.201101
## C1qa      2.2670625    0.972    0.396 6.629506e-13 7.937402e-15       1.509209
## C1qc      1.6955969    0.981    0.375 1.033899e-08 1.942429e-15       1.617483
## Cebpb     0.1487639    0.185    0.125 1.000000e+00 3.726173e-15       1.545982
## Sdc4      0.3247144    0.306    0.146 1.000000e+00 5.195913e-15       1.213547
## C1qb      1.9546166    0.972    0.417 1.524222e-10 8.047736e-12       1.113736
##       C24_pct.1 C24_pct.2 C24_p_val_adj
## Cd14      1.000     0.541  2.026193e-17
## C1qa      0.992     0.405  1.134493e-10
## C1qc      0.975     0.365  2.776313e-11
## Cebpb     0.983     0.581  5.325819e-11
## Sdc4      1.000     0.784  7.426519e-11
## C1qb      1.000     0.486  1.150263e-07
table(macrophages.3.integrated$sample.id)
## 
##  C0 C24 
## 156 192
colnames(conserved.markers.pos)
##  [1] "cluster"        "gene"           "max_pval"       "minimump_p_val"
##  [5] "marker.type"    "C0_p_val"       "C0_avg_log2FC"  "C0_pct.1"      
##  [9] "C0_pct.2"       "C0_p_val_adj"   "C24_p_val"      "C24_avg_log2FC"
## [13] "C24_pct.1"      "C24_pct.2"      "C24_p_val_adj"
library(dplyr)
top1.pos <- conserved.markers.pos %>% group_by(cluster) %>% top_n(n = 1, 
    wt = -minimump_p_val)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = -max_pval)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = C24_avg_log2FC)
head(top1.pos)
## # A tibble: 3 x 15
## # Groups:   cluster [3]
##   cluster gene  max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
##   <chr>   <chr>    <dbl>          <dbl> <chr>          <dbl>         <dbl>
## 1 Macro1  Cd14  3.90e-10       2.84e-21 POS         3.90e-10        1.99  
## 2 Macro2  Pfkp  1.52e- 2       2.28e-23 POS         1.52e- 2        0.0989
## 3 Macro3  Kap   4.60e- 6       7.21e-11 POS         4.60e- 6        1.77  
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## #   C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## #   C24_p_val_adj <dbl>
DefaultAssay(macrophages.3.integrated) <- "RNA"
f <- FeaturePlot(macrophages.3.integrated, features = top1.pos$gene, 
    min.cutoff = "q9", order = T, label = T, label.size = 8, 
    repel = T)
f

pdf(paste0(outdir, "/71_FeaturePlot_top_pos_markers_macrophages_3.pdf"), 
    width = 14, height = 14)
f
dev.off()
## png 
##   2

DotPlot of top markers

DefaultAssay(macrophages.3.integrated) <- "RNA"
d <- DotPlot(object = macrophages.3.integrated, features = top1.pos$gene, 
    cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.3.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/72_DotPlot_top_cell_types_pos_markers_clusters.macrophages_3.pdf"), 
    width = 7, height = 7)
d
dev.off()
## png 
##   2

DotPlot of top 10 markers

library(dplyr)
top10.pos <- conserved.markers.pos %>% group_by(cluster) %>% 
    top_n(n = 10, wt = -minimump_p_val)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10, 
    wt = -max_pval)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10, 
    wt = C24_avg_log2FC)
head(top10.pos)
## # A tibble: 6 x 15
## # Groups:   cluster [1]
##   cluster gene  max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
##   <chr>   <chr>    <dbl>          <dbl> <chr>          <dbl>         <dbl>
## 1 Macro1  Cd14  3.90e-10       2.84e-21 POS         3.90e-10         1.99 
## 2 Macro1  C1qa  7.94e-15       9.28e-17 POS         4.64e-17         2.27 
## 3 Macro1  C1qc  7.23e-13       3.88e-15 POS         7.23e-13         1.70 
## 4 Macro1  Cebpb 3.64e- 1       7.45e-15 POS         3.64e- 1         0.149
## 5 Macro1  Sdc4  4.12e- 2       1.04e-14 POS         4.12e- 2         0.325
## 6 Macro1  C1qb  8.05e-12       2.13e-14 POS         1.07e-14         1.95 
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## #   C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## #   C24_p_val_adj <dbl>
DefaultAssay(macrophages.3.integrated) <- "RNA"
d <- DotPlot(object = macrophages.3.integrated, features = top10.pos$gene, 
    cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.3.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/73_DotPlot_top10_cell_types_pos_markers_clusters.macrophages.3.pdf"), 
    width = 14, height = 7)
d
dev.off()
## png 
##   2

Segregation of clusters by various sources of uninteresting variation

Explore additional metrics, such as the number of UMIs and genes per cell, and mitochondrial gene expression by UMAP.

# Determine metrics to plot present in
# tcells.integrated@meta.data
metrics <- c("nCount_RNA", "nFeature_RNA", "percent.mito")

f <- FeaturePlot(macrophages.3.integrated, reduction = "umap", 
    features = metrics, pt.size = 0.4, order = TRUE, min.cutoff = "q10", 
    label = TRUE)
f

pdf(paste0(outdir, "/74_FeaturePlot_macrophages.3.integrated_umap_metrics.pdf"), 
    width = 14, height = 14)
f
dev.off()
## png 
##   2

Heatmap

macrophages.3.integrated <- ScaleData(macrophages.3.integrated, 
    verbose = T, features = rownames(macrophages.3.integrated))
d <- DoHeatmap(macrophages.3.integrated, features = top10.pos$gene, 
    group.by = "renamed.macrophages.3.subclusters", group.colors = colors.macrophages.3.subclusters) + 
    NoLegend()
d

pdf(paste0(outdir, "/75_DoHeatmap_top10_pos_macrophages_3.pdf"), 
    width = 15, height = 13)
d
dev.off()
## png 
##   2

DotPlot with known/useful markers

immune.cells <- c("Ptprc", "Spi1", "Cd3g")
b.cells <- c("Cd19", "Cd79a", "Cd79b")
t.cells <- "Cd3e"
nk.cells <- c("Ncr1")
myeloid.cells <- "Cd68"
neutrophils <- c("Csf3r", "S100a8", "S100a9", "Retnlg")
monocytes.macrophages.dendritic.cells <- "Csf1r"
dc <- c("Siglech", "Cd209a")
dc1 <- "Xcr1"
dc2 <- c("Ccl18", "Nr4a3", "Tnfsf12")
dc3 <- "Fscn1"
pdc <- "Tcf4"
monocytes <- c("Ly6c1", "Ly6c2", "Ccr2", "Cx3cr1", "Fcgr4")
macrophages <- c("Mrc1", "Lyz2", "Adgre1", "Axl", "Cxcl2", "Trem2", 
    "C1qc", "Fcgr1", "C5ar1")
am.markers <- c("Pparg", "Itgax", "Ear1", "Ear2", "Ear3", "Car4", 
    "Siglec5", "SiglecF")
im.markers <- c("Itgam", "Cx3cr1", "Ccl12", "Bhlme40", "Trem2", 
    "Cxcl12", "Csf1")
peribronchial.macrophages <- c("Ccr2", "Bhlhe40", "Bhlhe41")
perivascular.macrophages <- c("Lyve1", "F13a1")
proliferating.cells <- "Mki67"
mesenchymal.cells <- c("Col1a1", "Col1a2", "Pdgfra")
epithelial.cells <- c("Epcam", "Cdh1")
mast.cells <- "Mcpt8"
endothelial.cells <- c("Pecam1", "Emcn")
markers <- unique(c(immune.cells, b.cells, t.cells, nk.cells, 
    myeloid.cells, neutrophils, monocytes.macrophages.dendritic.cells, 
    dc, dc1, dc2, dc3, pdc, monocytes, macrophages, am.markers, 
    im.markers, peribronchial.macrophages, perivascular.macrophages, 
    proliferating.cells, mesenchymal.cells, epithelial.cells, 
    mast.cells, endothelial.cells))

Dot plot

DefaultAssay(macrophages.3.integrated) <- "RNA"
d <- DotPlot(object = macrophages.3.integrated, features = markers, 
    cols = c("#03539C", "#E82564"), dot.scale = 15, group.by = "renamed.macrophages.3.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/76_DotPlot_known_markers_macrophages_3.pdf"), 
    width = 18, height = 4)
d
dev.off()
## png 
##   2

Expression of genes of interest: Socs3, Il1r2, Il1rn

goi <- c("Socs3", "Il1r2", "Il1rn")
DefaultAssay(macrophages.3.integrated) <- "RNA"
f <- FeaturePlot(macrophages.3.integrated, features = goi, min.cutoff = "q9", 
    order = T, label = T, label.size = 8, repel = T)
f
pdf(paste0(outdir, "/77_FeaturePlot_goi_macrophages_3.pdf"), 
    width = 14, height = 14)
f
dev.off()
DefaultAssay(macrophages.3.integrated) <- "RNA"
d <- DotPlot(object = macrophages.3.integrated, features = goi, 
    cols = c("#03539C", "#E82564"), dot.scale = 15, group.by = "renamed.macrophages.3.subclusters") + 
    RotatedAxis()
d

pdf(paste0(outdir, "/78_DotPlot_goi_macrophages.pdf"), width = 6, 
    height = 4)
d
dev.off()
## png 
##   2

DIFFERENTIAL EXPRESSION ANALYSES

differential expression studies for the following groups: 1. C24 vs C0 within each cell type (from manual annotation annotation.1)

1. C24 vs C0 within each cell type (from manual annotation annotation.1)

macrophages.3.integrated$renamed.macrophages.3.subclusters.sample.id <- paste(macrophages.3.integrated$renamed.macrophages.3.subclusters, 
    macrophages.3.integrated$sample.id, sep = "_")
levels(factor(macrophages.3.integrated$renamed.macrophages.3.subclusters.sample.id))
## [1] "Macro1_C0"  "Macro1_C24" "Macro2_C0"  "Macro2_C24" "Macro3_C0" 
## [6] "Macro3_C24"
DefaultAssay(macrophages.3.integrated) <- "RNA"
output <- paste0(outdir, "/79_DE_macrophages_3_subclusters_C24_vs_C0_")
levels(factor(macrophages.3.integrated$renamed.macrophages.3.subclusters.sample.id))
## [1] "Macro1_C0"  "Macro1_C24" "Macro2_C0"  "Macro2_C24" "Macro3_C0" 
## [6] "Macro3_C24"
for (cell.type in levels(factor(macrophages.3.integrated$renamed.macrophages.3.subclusters))) {
    try({
        print(cell.type)
        ident1 <- paste0(cell.type, "_C24")
        ident2 <- paste0(cell.type, "_C0")
        Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
        condition.diffgenes <- FindMarkers(macrophages.3.integrated, 
            ident.1 = ident1, ident.2 = ident2, logfc.threshold = -Inf, 
            test.use = "wilcox", min.pct = -Inf, verbose = T, 
            min.cells.feature = -Inf, min.cells.group = -Inf)
        
        condition.diffgenes.adjpval.0.05 <- condition.diffgenes[condition.diffgenes$p_val_adj < 
            0.05, ]
        condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05[condition.diffgenes.adjpval.0.05$avg_log2FC > 
            0, ]
        condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05.upregulated.sorted[order(condition.diffgenes.adjpval.0.05.upregulated.sorted$avg_log2FC, 
            decreasing = T), ]
        condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05[condition.diffgenes.adjpval.0.05$avg_log2FC < 
            0, ]
        condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05.downregulated.sorted[order(condition.diffgenes.adjpval.0.05.downregulated.sorted$avg_log2FC, 
            decreasing = F), ]
        list_of_datasets <- list(sig.upregulated.genes.ordered.by.avg.logFC = condition.diffgenes.adjpval.0.05.upregulated.sorted, 
            sig.downregulated.genes.ordered.by.avg.logFC = condition.diffgenes.adjpval.0.05.downregulated.sorted, 
            all.genes.ordered.by.adj.pval = condition.diffgenes)
        openxlsx::write.xlsx(list_of_datasets, file = paste0(output, 
            cell.type, ".xlsx"), row.names = T)
    })
}
## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"
for (cell.type in levels(factor(macrophages.3.integrated$renamed.macrophages.3.subclusters))) {
    try({
        print(cell.type)
        filename <- paste0(output, cell.type, ".xlsx")
        sheets <- openxlsx::getSheetNames(filename)
        SheetList <- lapply(sheets, openxlsx::read.xlsx, xlsxFile = filename)
        names(SheetList) <- sheets
        
        condition.diffgenes.adjpval.0.05.upregulated.sorted <- SheetList[[1]]
        rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted) <- condition.diffgenes.adjpval.0.05.upregulated.sorted[, 
            1]
        condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05.upregulated.sorted[, 
            -1]
        
        condition.diffgenes.adjpval.0.05.downregulated.sorted <- SheetList[[2]]
        rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted) <- condition.diffgenes.adjpval.0.05.downregulated.sorted[, 
            1]
        condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05.downregulated.sorted[, 
            -1]
        
        
        Idents(annotated) <- "renamed.macrophages.3.subclusters"
        
        try({
            f <- FeaturePlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1], 
                split.by = "sample.id", cols = c("grey", "#E82564"), 
                order = T, min.cutoff = "q10", max.cutoff = "q90")
            pdf(paste0(output, cell.type, "_FeaturePlot_top_upregulated_split_by_sample.pdf"), 
                width = 14, height = 7)
            print(f)
            dev.off()
        })
        
        try({
            f <- FeaturePlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1], 
                split.by = "sample.id", cols = c("grey", "#E82564"), 
                order = T, min.cutoff = "q10", max.cutoff = "q90")
            pdf(paste0(output, cell.type, "_FeaturePlot_top_downregulated_split_by_sample.pdf"), 
                width = 14, height = 7)
            print(f)
            dev.off()
        })
        
        try({
            plot <- VlnPlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1], 
                group.by = "renamed.macrophages.3.subclusters", 
                split.by = "sample.id", pt.size = 0, combine = FALSE, 
                cols = colors.samples)
            pdf(paste0(output, cell.type, "_VlnPlot_top_upregulated_group_by_cluster_split_by_sample.pdf"), 
                width = 11, height = 8.5)
            print(plot)
            dev.off()
        })
        
        try({
            plot <- VlnPlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1], 
                group.by = "renamed.macrophages.3.subclusters", 
                split.by = "sample.id", pt.size = 0, combine = FALSE, 
                cols = colors.samples)
            pdf(paste0(output, cell.type, "_VlnPlot_top_downregulated_group_by_cluster_split_by_sample.pdf"), 
                width = 11, height = 8.5)
            print(plot)
            dev.off()
        })
        
        ## DotPlot with 50 most upregulated genes
        Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
        try({
            d <- DotPlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1:50], 
                cols = c("#03539C", "#E82564"), dot.scale = 8, 
                group.by = "renamed.macrophages.3.subclusters.sample.id") + 
                RotatedAxis()
            pdf(paste0(output, cell.type, "_DotPlot_top50_upregulated_groupby_annotation.1.sample.id.pdf"), 
                width = (4 + min(50, length(rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted))) * 
                  15/50), height = 17)
            print(d)
            dev.off()
        })
        
        ## DotPlot with 50 most downregulated genes
        Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
        try({
            d <- DotPlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1:50], 
                cols = c("#03539C", "#E82564"), dot.scale = 8, 
                group.by = "renamed.macrophages.3.subclusters.sample.id") + 
                RotatedAxis()
            pdf(paste0(output, cell.type, "_DotPlot_top50_downregulated_groupby_annotation.1.sample.id.pdf"), 
                width = (4 + min(50, length(rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted))) * 
                  15/50), height = 17)
            print(d)
            dev.off()
        })
    })
}
## [1] "Macro1"
## [[1]]
## 
## [[1]]
## [1] "Macro2"
## [[1]]
## 
## [[1]]
## [1] "Macro3"
## Error : None of the requested features were found: NA in slot data
## [[1]]
## 
## Error in FetchData(object = object, vars = features, slot = slot) : 
##   None of the requested variables were found: NA
## Error in h(simpleError(msg, call)) : 
##   error in evaluating the argument 'x' in selecting a method for function 'mean': non-numeric argument to mathematical function

Plots for genes of interest

for (gene in goi) {
    try({
        Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
        try({
            f <- FeaturePlot(macrophages.3.integrated, features = gene, 
                split.by = "sample.id", cols = c("grey", "#E82564"), 
                order = T, min.cutoff = "q10", max.cutoff = "q90", 
                label = T)
            pdf(paste0(outdir, "/80_", gene, "_FeaturePlot_split_by_sample.pdf"), 
                width = 14, height = 7)
            print(f)
            dev.off()
        })
        
        try({
            plot <- VlnPlot(macrophages.3.integrated, features = gene, 
                group.by = "renamed.macrophages.3.subclusters", 
                split.by = "sample.id", pt.size = 0, combine = FALSE, 
                cols = colors.samples)
            pdf(paste0(outdir, "/81_", gene, "_VlnPlot_group_by_cluster_split_by_sample.pdf"), 
                width = 11, height = 8.5)
            print(plot)
            dev.off()
        })
    })
}
## [[1]]
## [[1]]
## [[1]]
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
try({
    d <- DotPlot(macrophages.3.integrated, features = goi, cols = c("#03539C", 
        "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.3.subclusters.sample.id") + 
        RotatedAxis()
    pdf(paste0(outdir, "/82_", paste(goi, sep = "_"), "_DotPlot_groupby_annotation.1.sample.id.pdf"), 
        width = (4 + length(goi) * 15/50), height = 17)
    print(d)
    dev.off()
})
## png 
##   2

Plots for more genes of interest: F480 (antibody), CD11b, CD11c, CD64, MHCII

In mouse: F4/80 (antibody) > Emr1, Adgre1 CD11b -> Itgam CD11c -> Itgax CD64 -> Fcgr1 MHCII -> H2-Ab

goi <- c("Adgre1", "Itgam", "Itgax", "Fcgr1", "H2-Ab1")
DefaultAssay(macrophages.3.integrated) <- "RNA"
for (gene in goi) {
    c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.", 
        node_colour = gene, node_colour_aggr = "mean")
    pdf(paste0(outdir, "/83_clustree_", gene, "_macrophages.pdf"))
    print(c)
    dev.off()
}
for (gene in goi) {
    try({
        Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
        try({
            f <- FeaturePlot(macrophages.3.integrated, features = gene, 
                split.by = "sample.id", cols = c("grey", "#E82564"), 
                order = T, min.cutoff = "q10", max.cutoff = "q90", 
                label = T)
            pdf(paste0(outdir, "/84_", gene, "_FeaturePlot_split_by_sample.pdf"), 
                width = 14, height = 7)
            print(f)
            dev.off()
        })
        
        try({
            plot <- VlnPlot(macrophages.3.integrated, features = gene, 
                group.by = "renamed.macrophages.3.subclusters", 
                split.by = "sample.id", pt.size = 0, combine = FALSE, 
                cols = colors.samples)
            pdf(paste0(outdir, "/85_", gene, "_VlnPlot_group_by_cluster_split_by_sample.pdf"), 
                width = 11, height = 8.5)
            print(plot)
            dev.off()
        })
    })
}
## [[1]]
## [[1]]
## [[1]]
## [[1]]
## [[1]]
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
try({
    d <- DotPlot(macrophages.3.integrated, features = goi, cols = c("#03539C", 
        "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.3.subclusters.sample.id") + 
        RotatedAxis()
    pdf(paste0(outdir, "/86_", paste(goi, sep = "_"), "_DotPlot_groupby_macrophages_subclusters.sample.id.pdf"), 
        width = (4 + length(goi) * 15/50), height = 17)
    print(d)
    dev.off()
})
## png 
##   2

Session Information

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.8 (Santiago)
## 
## Matrix products: default
## BLAS:   /gpfs/fs1/data/omicscore/Privratsky-Privratsky-20210215/scripts/conda/envs/privratsky/lib/libblas.so.3.8.0
## LAPACK: /gpfs/fs1/data/omicscore/Privratsky-Privratsky-20210215/scripts/conda/envs/privratsky/lib/liblapack.so.3.8.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] Polychrome_1.2.6            clustree_0.4.3             
##  [3] ggraph_2.0.4                scater_1.18.0              
##  [5] SingleCellExperiment_1.12.0 data.table_1.13.6          
##  [7] gridExtra_2.3               forcats_0.5.1              
##  [9] stringr_1.4.0               purrr_0.3.4                
## [11] readr_1.4.0                 tidyr_1.1.2                
## [13] tibble_3.0.6                tidyverse_1.3.0            
## [15] SingleR_1.4.0               celldex_1.0.0              
## [17] SummarizedExperiment_1.20.0 Biobase_2.50.0             
## [19] GenomicRanges_1.42.0        GenomeInfoDb_1.26.0        
## [21] IRanges_2.24.0              S4Vectors_0.28.0           
## [23] BiocGenerics_0.36.0         MatrixGenerics_1.2.0       
## [25] matrixStats_0.58.0          ggsci_2.9                  
## [27] cowplot_1.1.1               plyr_1.8.6                 
## [29] patchwork_1.1.1             dplyr_1.0.4                
## [31] ggplot2_3.3.3               SeuratObject_4.0.0         
## [33] Seurat_4.0.0               
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.1.4                    reticulate_1.18              
##   [3] tidyselect_1.1.0              RSQLite_2.2.4                
##   [5] AnnotationDbi_1.52.0          htmlwidgets_1.5.3            
##   [7] grid_4.0.3                    BiocParallel_1.24.0          
##   [9] Rtsne_0.15                    munsell_0.5.0                
##  [11] codetools_0.2-18              ica_1.0-2                    
##  [13] future_1.21.0                 miniUI_0.1.1.1               
##  [15] withr_2.4.1                   colorspace_2.0-0             
##  [17] highr_0.8                     knitr_1.31                   
##  [19] rstudioapi_0.13               ROCR_1.0-11                  
##  [21] tensor_1.5                    gbRd_0.4-11                  
##  [23] listenv_0.8.0                 Rdpack_2.1                   
##  [25] labeling_0.4.2                GenomeInfoDbData_1.2.4       
##  [27] polyclip_1.10-0               farver_2.0.3                 
##  [29] bit64_4.0.5                   parallelly_1.23.0            
##  [31] vctrs_0.3.6                   generics_0.1.0               
##  [33] xfun_0.20                     BiocFileCache_1.14.0         
##  [35] R6_2.5.0                      graphlayouts_0.7.1           
##  [37] ggbeeswarm_0.6.0              rsvd_1.0.3                   
##  [39] bitops_1.0-6                  spatstat.utils_2.0-0         
##  [41] cachem_1.0.4                  DelayedArray_0.16.0          
##  [43] assertthat_0.2.1              promises_1.2.0.1             
##  [45] scales_1.1.1                  beeswarm_0.2.3               
##  [47] gtable_0.3.0                  beachmat_2.6.0               
##  [49] globals_0.14.0                goftest_1.2-2                
##  [51] tidygraph_1.2.0               rlang_0.4.10                 
##  [53] scatterplot3d_0.3-41          splines_4.0.3                
##  [55] lazyeval_0.2.2                checkmate_2.0.0              
##  [57] broom_0.7.5                   BiocManager_1.30.10          
##  [59] yaml_2.2.1                    reshape2_1.4.4               
##  [61] abind_1.4-5                   modelr_0.1.8                 
##  [63] backports_1.2.1               httpuv_1.5.5                 
##  [65] tools_4.0.3                   ellipsis_0.3.1               
##  [67] RColorBrewer_1.1-2            ggridges_0.5.3               
##  [69] Rcpp_1.0.6                    sparseMatrixStats_1.2.0      
##  [71] zlibbioc_1.36.0               RCurl_1.98-1.2               
##  [73] ps_1.5.0                      rpart_4.1-15                 
##  [75] deldir_0.2-9                  viridis_0.5.1                
##  [77] pbapply_1.4-3                 zoo_1.8-8                    
##  [79] haven_2.3.1                   ggrepel_0.9.1                
##  [81] cluster_2.1.1                 fs_1.5.0                     
##  [83] magrittr_2.0.1                RSpectra_0.16-0              
##  [85] scattermore_0.7               openxlsx_4.2.3               
##  [87] lmtest_0.9-38                 reprex_1.0.0                 
##  [89] RANN_2.6.1                    fitdistrplus_1.1-3           
##  [91] hms_1.0.0                     mime_0.10                    
##  [93] evaluate_0.14                 xtable_1.8-4                 
##  [95] readxl_1.3.1                  compiler_4.0.3               
##  [97] KernSmooth_2.23-18            crayon_1.4.1                 
##  [99] htmltools_0.5.1.1             mgcv_1.8-33                  
## [101] later_1.1.0.1                 lubridate_1.7.10             
## [103] DBI_1.1.1                     formatR_1.7                  
## [105] tweenr_1.0.1                  ExperimentHub_1.16.0         
## [107] dbplyr_2.1.0                  MASS_7.3-53.1                
## [109] rappdirs_0.3.3                Matrix_1.3-2                 
## [111] cli_2.3.0                     rbibutils_2.0                
## [113] metap_1.1                     igraph_1.2.6                 
## [115] pkgconfig_2.0.3               plotly_4.9.3                 
## [117] scuttle_1.0.0                 xml2_1.3.2                   
## [119] vipor_0.4.5                   XVector_0.30.0               
## [121] rvest_1.0.0                   digest_0.6.27                
## [123] sctransform_0.3.2             RcppAnnoy_0.0.18             
## [125] spatstat.data_2.0-0           rmarkdown_2.6                
## [127] cellranger_1.1.0              leiden_0.3.7                 
## [129] uwot_0.1.10                   DelayedMatrixStats_1.12.0    
## [131] curl_4.3                      shiny_1.6.0                  
## [133] lifecycle_1.0.0               nlme_3.1-152                 
## [135] jsonlite_1.7.2                BiocNeighbors_1.8.0          
## [137] fansi_0.4.2                   limma_3.46.0                 
## [139] viridisLite_0.3.0             pillar_1.4.7                 
## [141] lattice_0.20-41               fastmap_1.1.0                
## [143] httr_1.4.2                    survival_3.2-7               
## [145] interactiveDisplayBase_1.28.0 glue_1.4.2                   
## [147] zip_2.1.1                     spatstat_1.64-1              
## [149] png_0.1-7                     BiocVersion_3.12.0           
## [151] bit_4.0.4                     ggforce_0.3.2                
## [153] stringi_1.5.3                 blob_1.2.1                   
## [155] BiocSingular_1.6.0            AnnotationHub_2.22.0         
## [157] memoise_2.0.0                 irlba_2.3.3                  
## [159] future.apply_1.7.0
writeLines(capture.output(sessionInfo()), "./scripts/9_Subcluster_Macrophages/9_Subcluster_Macrophages.sessionInfo.txt")